! Questions for Laurent
! 1. IOQFieldValueVar. Try to explain why need this when are taxes on outputs 
!      but not when are taxes on inputs.
!   This must be documented. I think that
!    (a) in o: field, the q: field coeff is a pre-sim value INCLUDING taxes but
!    (b) in i: field, the q: field coeff is a pre-sim value NOT INCLUDING taxes.
!    Check this with Laurent. Put some of this into the labelinfo.
! 2. Need to test this with 2 taxes on i: line. In particular need to 
!     test mixedTAB output in that case.
! 3. Need an example with just one tax on o: line. This is needed to test correct
!     use of IOQFieldValueVar in that case.
! 4. $demand
!   (a) Always exactly one d: field? [Check with Laurent]
!   (b) What happens if are arguments involved in d: fields?
! 5. Tax revenue. Can a: field in i: or o: field have arguments? I assume so. example??
!     In adding revenue from a: terms, I have assumed no arguments. Will probably
!      need to change this.
! Documentation.
! 1. Rule for read header:
!   (a) if length <= 4, try to use it
!   (b) otherwise omit vowels from right-hand end (though not the first letter if it is a vowel)
!        to get a candidate. Stop omitting vowels if length gets to 4.
!   (c) Truncate candidate from (a) or (b) to length 4.
!   (d) Use this is not already used.
!   (e) Otherwise ?? (nothing implemented yet)
! 2. Check all calls to ErrorStop for user errors
! 3. Check all calls to GiveWarning for user advice  

!

module mpsGP

! module used in conversion to GEMPACK TAB and Command files

	use mpsge

! Whether to translate to a mixed or a linear TAB file
!  Only one of these will be true
	logical :: MixedTABFile, LinearTABFile

! Lengths
	integer, parameter :: llst=12  ! Max length of set name
	integer, parameter :: llcs=12  ! Max length of Coefficient name
	integer, parameter :: llvc=15  ! Max length of (linear) Variable name
	integer, parameter :: llrc=15  ! Max length of real constant
	integer, parameter :: llindex=12  ! Max length of index
	
! Declarations/Types added by KP

	integer :: NumSets  !! Number of sets in the MPS file

	type SetGP
	  character (len=llst) :: name
	  type(SetGP), pointer :: next  !! pointer to the next setGP in the list
	end type SetGP

	type(setGP), pointer :: SetGPList !! List of all the SetGPs
	type(setGP), pointer :: SetGPLast !! Last of all the SetGPs

! For subset mappings
	logical, dimension(:,:), allocatable :: IsSubsetGP  !! Tells if each set is a subset of any other set

	type GPCoeffVarReal    ! GEMPACK Coefficient or Variable or Real
	  character :: type  ! 'c' for coefficient or 'v' for variable or 'r' for real
	  character (len=MAX(llcs,llvc,llrc)) :: name    ! Coeff/Var name or Real as string
	  type (domain), pointer :: arg      ! Sets over which arguments range
	  character (len=80) :: labelinfo  ! labelling information (between #)
! Next two only apply to a variable
	  character :: VarType  ! 'p' for %-change, 'c' for change
	  character (len=12) :: OrigLevel  ! Coeff or Real which is the ORIG_LEVEL ! ignore if OrigLevel is ' '
! Next only applies to a coefficient
	  logical :: IsParam  ! Tells if this coefficient is a parameter
	end type GPCoeffVarReal

        integer, parameter :: iotabgp=11   !! Unit number of the output GEMPACK TAB file
        integer, parameter :: iocmfgp=12   !! Unit number of the output GEMPACK Command file
        character (len=255) :: tabfilename, cmffilename, tabfileroot

	integer, parameter :: llhd=4   !! Length of a header in GEMPACK
	type HeaderGP  ! List of headers
	  character (len=llhd) :: header
	  type (HeaderGP), pointer :: next  !! Pointer to the next one in the list
	end type HeaderGP

	type (HeaderGP), pointer :: ReadHeaderGPList !! List of headers used for reading
	type (HeaderGP), pointer :: ReadHeaderGPLast !! Last header in list of headers used for reading

	integer LengthGPDeclarationName
	parameter (LengthGPDeclarationName = 20)
	type GPDeclaration   ! List of various types of GEMPACK declarations
	  character (len=2) :: type  !! 'cs' for Coeff, 'vc' for Variable, 'eq' for equation
	  character (len=LengthGPDeclarationName) :: name  !! Coeff/Var/Equation name
	  type (GPDeclaration), pointer :: next !! Pointer to the next declaration in this list
	end type GPDeclaration

	type (GPDeclaration), pointer :: GPDeclarationList !! Pointer to first in GPDeclarationList 
	type (GPDeclaration), pointer :: GPDeclarationLast !! Pointer to last in GPDeclarationList

	type GPName   ! List of various types of GEMPACK declarations
	  character (len=LengthGPDeclarationName) :: name  !! Coeff/Var/Equation name
	  character :: UserSystemType !! 'u' for names supplied by user, 's' for names made up by our system
	  type (domain), pointer :: arg      ! Sets over which arguments range
! next only for variables
	  character :: type ! 'v' for variable, 'c' for coefficient
	  character :: ExogEndogStatus ! 'n' for endogenous, 'x' for exogenous, '?' for not sure
	  type (GPName), pointer :: next !! Pointer to the next name in this list
	end type GPName

	type (GPName), pointer :: GPNameList !! Pointer to first in GPNameList 
	type (GPName), pointer :: GPNameLast !! Pointer to last in GPNameList

	type GPRead   ! List of names of Coefficients whose values have been read
	  character (len=LengthGPDeclarationName) :: name  !! Coeff name
! ??	  type (domain), pointer :: arg      ! Sets over which arguments range
	  type (GPRead), pointer :: next !! Pointer to the next name in this list
	end type GPRead

	type (GPRead), pointer :: GPReadList !! Pointer to first in GPReadList 
	type (GPRead), pointer :: GPReadLast !! Pointer to last in GPReadList

	type (GPCoeffVarReal) :: WalrasSlackVar !! Variable determined by difference between 2 sides of omitted (Walras) equation. Should be zero.

        type functionGP
! Additions to Tom Rutherford's type function. These additions needed for GEMPACK
! This is actually a list of functionGP's
	  type (function), pointer :: functionMPSGE  !! Pointer to the MPSGE information

	  logical :: IsProdFunctionGP, IsDemandFunctionGP

!  The first lot below only apply for a $prod function
	  type (GPCoeffVarReal) :: ActivityLevelVar !! Variable after %prod:
	  logical :: SElastPresent, TElastPresent, SElastIsConstant, TElastIsConstant
	  type (field), pointer :: SElastField, TElastField
	  type (GPCoeffVarReal) :: SElasticity, TElasticity

! Next are for $prod functions
	  type (GPCoeffVarReal) :: TotalRevCoeff !! Total revenue as coefficient (relates to outputs)
	  type (GPCoeffVarReal) :: MargRevVar !! Marginal revenue as variable (relates to outputs)
	  type (GPCoeffVarReal) :: TotalCostCoeff !! Total cost as coefficient (relates to inputs)
	  type (GPCoeffVarReal) :: MargCostVar !! Marginal cost as variable (relates to inputs)
	  integer :: NumIFields, NumOFields  !! Number of i: and o: fields

	  type (netputGP), pointer :: FirstNetputGP   !! Points to first netputGP in netputGPList associated with this function 

! Next are for $demand functions
	  integer :: NumDFields, NumEFields  !! Number of d: and e: fields
	  type (GPCoeffVarReal) :: IncomeVar !! %-change in value of income for this agent
	  type (GPCoeffVarReal) :: TotalIncomeCoeff !! Total income for this consumer
		! Is equal to total value of endowments (summed over all e: lines) plus total tax income

	  type (endowmentGP), pointer :: FirstEndowmentGP   !! Points to first endowmentGP in endowmentGPList associated with this function 

          type (functionGP), pointer :: next      !! Pointer in linked list of all functionGPs
        end type functionGP

        type (functionGP), pointer :: FunctionGPList !! Pointer to the first functionGP in linked list
        type (functionGP), pointer :: FunctionGPLast !! Pointer to the last functionGP in linked list
	
        type netputGP
! Additions to Tom Rutherford's type netput. These additions needed for GEMPACK
! This is actually a list of netputGP's
	  type(netput), pointer :: netputMPSGE !! Pointer to the MPSGE netput

! The GEMPACK arguments for this netputGP may be different from netputMPSGE%arg. The
!  sets will be the same but they may be in the same order. For example,
! $prod:y(s)          t:etrn(s)
!  o:p(o)            q:supply(s,o)       p:0.8  a:gov  t:0.2     a:gov  t:t_p(o)
!  i:p(o)            q:interm(o,s)
!  i:pva(s)          q:totfac
! For netput beginning o:p(o), the netputMPSGE%arg are (s,o). Since the arguments in the q: field
!   are in the same order, IOArgGP will be the same (s,o).
! But for netput beginning i:p(o), the netputMPSGE%arg is (s,o). However, since the arguments in
!  the q: field are there sets but in a different order (ie, the q: field parameter has "full rank"),
!  IOArgGP will be (o,s) in that order.
! At present, IOArgGP and IOQFieldValueCoeff%arg are always the same. But use IOArgGP in the code
!  in case this changes in some special cases which are excluded at present.
	  type (domain), pointer :: IOArgGP  

!		IOPriceVar is the market price
	  type (GPCoeffVarReal) :: IOPriceVar    !! %-change in price (The price after o: or i: at start of line) 
	  type (GPCoeffVarReal) :: IOPriceCoeff    !! price (levels value) (The price after o: or i: at start of line) 
	  type (GPCoeffVarReal) :: IOQFieldQuantityCoeff  !! quantity (level) of each input or output (q: field)  
	  type (GPCoeffVarReal) :: IOQFieldQuantityVar  !! %-change in input or output quantity
	  type (GPCoeffVarReal) :: IOQFieldValueCoeff  !! value (level) of input or output
	  type (GPCoeffVarReal) :: IOQFieldValueVar    !! %-change in value of input or output
! Names given for this example
! $prod:y(s)          t:etrn(s)
!  o:p(o)            q:200
! IOPriceVar is p(o) and IOPriceCoeff is P_L(o)
! IOQFieldQuantityCoeff is QS_Y_P(s,o) and IOQFieldQuantityVar is p_qs_y_p(s,o)
! IOQFieldValueCoeff is VS_Y_P(s,o) and IOQFieldValueVar is p_vs_y_p(s,o)


! Next will be a real (eg q:200) or a Coefficient (eg q:totfac)
! Note that IOQFieldQuantityCoeffOrRealMPSGE is not set up if QFieldIsExpression is true
	  type (GPCoeffVarReal) :: IOQFieldQuantityCoeffOrRealMPSGE  !! MPSGE value (initial level) of input or output
! Now handle case where is an experssion in the q: field
	  logical :: QFieldIsExpression  !! Tells if QField is an expression
	  character (len=512) :: QFieldExpressionString   !! The expression in the q: field (if is such) 

	  integer :: NumTaxTypes
! Next only needed if are one or more taxing agents
	  type (GPCoeffVarReal) :: IOPriceIncTaxesVar    !! %-change in price including taxes
! Next only needed for outputs when are one or more taxing agents. 
!		[Need this for update of the values in the output case. Need division in levels
!			equation for this value, so can't use product update. For inputs can use a product update.
!		        The difference relates to whether value includes taxes or not. Ken to check with Laurent to
!				try to express and explain this better.]

	  type (taxGP), pointer :: FirstTaxGP   !! Points to first taxGP in TaxGPList associated with this netput 

! Next only needed if are 2 or more tax types for this netput
	  type (GPCoeffVarReal) :: TotalTaxPowerCoeff  !! power of tax (all tax types)
	  type (GPCoeffVarReal) :: TotalTaxPowerVar  !! %-change in power of tax (all tax types)

	  type (netputGP), pointer :: next      !! Pointer in linked list of all netputGPs
        end type netputGP

        type (netputGP), pointer :: NetputGPList !! Pointer to the first netputGP in linked list
        type (netputGP), pointer :: NetputGPLast !! Pointer to the last netputGP in linked list

	type taxGP
! Additions to Tom Rutherford's type taxlist. These additions needed for GEMPACK
! This is actually a list of taxlistGP's
	  type(taxlist), pointer :: taxlistMPSGE !! Pointer to the MPSGE taxlist

	  type (GPCoeffVarReal) :: TaxPowerCoeff  !! power of tax
	  type (GPCoeffVarReal) :: TaxRateCoeff  !! rate of tax (number like 0.13)
	  type (GPCoeffVarReal) :: TaxRevenueCoeff  !! tax revenue
	  type (GPCoeffVarReal) :: TaxPowerVar  !! %-change in power of tax
	  type (GPCoeffVarReal) :: TaxRateChangeVar  !! change in tax rate (number like 0.05)
	  type (GPCoeffVarReal) :: TaxRevenueChangeVar  !! change in tax revenue
! Next will be a real (eg t:0.2) or a Coefficient (eg t:t_p(o))
	  type (GPCoeffVarReal) :: TaxRateCoeffOrRealMPSGE  !! MPSGE value (initial level) of tax rate

	  type(taxGP), pointer :: next    !! Pointer to next in linked list of all taxGPs
	end type taxGP

        type (taxGP), pointer :: TaxGPList !! Pointer to the first taxGP in linked list
        type (taxGP), pointer :: TaxGPLast !! Pointer to the last taxGP in linked list

	type endowmentGP
! Additions to Tom Rutherford's type endowment. These additions needed for GEMPACK
! This is actually a list of endowmentGP's
	  type(endowment), pointer :: endowmentMPSGE !! Pointer to the MPSGE endowment

! EndowArgGP is analagous to IOArgGP
! At present, EndowArgGP and EndowQFieldValueCoeff%arg are always the same. But use EndowArgGP in the code
!  in case this changes in some special cases which are excluded at present.
	  type (domain), pointer :: EndowArgGP  

! The next 6 are analagous to the corresponding IO things in a netputGP
!			EndowPriceVar is the market price of this endowment
	  type (GPCoeffVarReal) :: EndowPriceVar    !! %-change in price (The price after e: at start of line)
	  type (GPCoeffVarReal) :: EndowPriceCoeff    !! price (levels value) (The price after e: at start of line)
	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeff  !! quantity (level) of each endowment q: field  
	  type (GPCoeffVarReal) :: EndowQFieldQuantityVar  !! change in quantity of each endowment (q: field) 
	  type (GPCoeffVarReal) :: EndowQFieldValueCoeff  !! value (level) of endowment
	  type (GPCoeffVarReal) :: EndowQFieldValueVar  !! %-change in value (level) of endowment
! Names are given for the example:
! $demand:ra
! e:plab  q:200
! EndowPriceVar is plab and EndowPriceCoeff is PLAB_L
! EndowQFieldQuantityCoeff is QS_RA_PLAB and EndowQFieldQuantityVar is p_qs_ra_plab [or c_qs_ra_plab]
! EndowQFieldValueCoeff is VS_RA_PLAB and EndowQFieldValueVar is p_vs_ra_plab [or c_vs_ra_plab]

! Next will be a real (eg q:200) or a Coefficient (eg q:totfac)
! Note that EndowQFieldQuantityCoeffOrRealMPSGE is not set up if QFieldIsExpression is true
	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeffOrRealMPSGE  !! MPSGE value (initial level) of endowment

! Now handle case where is an experssion in the q: field
	  logical :: QFieldIsExpression  !! Tells if QField is an expression
	  character (len=512) :: QFieldExpressionString   !! The expression in the q: field (if is such) 

	  logical :: HasRField
! Next only set up if HasRField is true. The GEMPACK name for EndowRFieldCoeff is the same as in the MGE file.
! EndowRFieldVar is always a change variable with name c_<EndowRFieldCoeff>.
	  type (GPCoeffVarReal) :: EndowRFieldCoeff  !! level of r: variable
	  type (GPCoeffVarReal) :: EndowRFieldVar  !! change in level of r: variable

! EndowUnRationedQuantityCoeff/Var are of no economic interest - we have them purely to implement the MPSGE
!   idea of "rationing" 
	  type (GPCoeffVarReal) :: EndowUnrationedQuantityCoeff  !! unrationed quantity of endowment (level)
	  type (GPCoeffVarReal) :: EndowUnrationedQuantityVar  !! change in unrationed quantity of endowment
! Names are given for the example:
! $demand:ra
! e:plab  q:200  r:rls
! EndowPriceVar is plab
! EndowQFieldQuantityCoeff is QS_RA_PLAB and EndowQFieldQuantityVar is p_qs_ra_plab [or c_qs_ra_plab]
! EndowQFieldValueCoeff is VS_RA_PLAB and EndowQFieldValueVar is p_vs_ra_plab [or c_vs_ra_plab]
! EndowRFieldCoeff is RLS and EndowRFieldVar is c_RLS
! EndowUnrationedQuantityCoeff is US_RA_PLAB and EndowUnrationedQuantityVar is c_us_ra_plab.

! [04/04/04] These will be omitted
!	  type (GPCoeffVarReal) :: EndowRationedQuantityVar  !! %-change in quantity of each endowment  
!	  type (GPCoeffVarReal) :: EndowRationedQFieldValueCoeff  !! value (level) of endowment
!	  type (GPCoeffVarReal) :: EndowRationedQFieldValueVar  !! %-change in value (level) of endowment

	  type(endowmentGP), pointer :: next    !! Pointer to next in linked list of all endowmentGPs
	end type endowmentGP

        type (endowmentGP), pointer :: EndowmentGPList !! Pointer to the first endowmentGP in linked list
        type (endowmentGP), pointer :: EndowmentGPLast !! Pointer to the last endowmentGP in linked list

contains

	subroutine AddToGPNameList( Name, UserSystemType, CoeffVarType, Args)
! This adds at the end of GPNameList a GPName with the specified
!		Name, UserSystemType and Args
! Input
	character (len=*) Name, UserSystemType, CoeffVarType
	type (domain), pointer :: Args      ! Sets over which arguments range
! Local
	type (GPName), pointer :: ThisGPName

	allocate(ThisGPName)
	ThisGPName%name = Name
	nullify(ThisGPName%next)
! Complete info in ThisGPName
	ThisGPName%UserSystemType = UserSystemType
	ThisGPName%type = CoeffVarType
	ThisGPName%arg => Args
	if(CoeffVarType.eq.'v') then
! Assume originally that all are exogenous. Mark ones as endogenous when write equation for them
	  ThisGPName%ExogEndogStatus = 'x'
	else
	  ThisGPName%ExogEndogStatus = ' '
	end if

	if(.not.associated(GPNameList)) then
	  GPNameList => ThisGPName
	  GPNameLast => ThisGPName
	else
	  GPNameLast%next => ThisGPName
	  GPNameLast => ThisGPName
	end if
	nullify(GPNameLast%next)

	return
	end subroutine AddToGPNameList

	subroutine AddToReadHeaderGPList( ThisHeader)
! This adds at the end of ReadHeaderGPListList a header
! Input
	character (len=*) ThisHeader
! Local
	type (HeaderGP), pointer :: NewHeaderGP

	allocate(NewHeaderGP)
	NewHeaderGP%Header = ThisHeader
	nullify(NewHeaderGP%next)

	if(.not.associated(ReadHeaderGPList)) then
	  ReadHeaderGPList => NewHeaderGP
	  ReadHeaderGPLast => NewHeaderGP
	else
	  ReadHeaderGPLast%next => NewHeaderGP
	  ReadHeaderGPLast => NewHeaderGP
	end if
	nullify(ReadHeaderGPLast%next)

	return
	end subroutine AddToReadHeaderGPList

	logical function IsInReadHeaderGPList(ThisString)
! Input
	character (len=*), intent(in) :: ThisString
! Local
	type (HeaderGP), pointer :: CurrentHeaderGP
	character (len=llhd) :: AbbrevString, Head2

	AbbrevString = ThisString

	CurrentHeaderGP => ReadHeaderGPList
	do while(associated(CurrentHeaderGP))
	  Head2 = CurrentHeaderGP%Header
!	  write(*,*)'IsInReadHeaderGPList. AbbrevString='//trim(AbbrevString)// &
!		' Head2 =',trim(head2)
	  if(EqualStringsIgnoreCase(AbbrevString, Head2)) then
	    IsInReadHeaderGPList = .true.
	    Return
	  end if

	  CurrentHeaderGP => CurrentHeaderGP%next
	end do

	IsInReadHeaderGPList = .false.
	return
	end function IsInReadHeaderGPList

	integer function LastVowel(ThisString)
! Marks position of last vowel in ThisString
! Input
	character (len=*), intent(in) :: ThisString
! Local
	character (len=512) :: TempString
	integer n1,n2,n3,n4,n5

	TempString = ThisString
	call lower(TempString)
	n1=INDEX(TempString, 'a', back=.true.)
	n2=INDEX(TempString, 'e', back=.true.)
	n3=INDEX(TempString, 'i', back=.true.)
	n4=INDEX(TempString, 'o', back=.true.)
	n5=INDEX(TempString, 'u', back=.true.)
	LastVowel = MAX(n1,n2,n3,n4,n5)
!	write(*,*)'LastVowel. ThisString=', trim(ThisString), ' Result=', LastVowel
	
	return
	end function LastVowel

	subroutine GetReadHeader(NameToRead, NewHeader)
! Input
	character (len=*), intent(in) :: NameToRead
! Output
	character (len=*), intent(out) :: NewHeader
! Must avoid any headers in ReadHeaderGPList
! Local
	type (HeaderGP), pointer :: CurrentHeader
	character (len=llhd) :: TryHead1, TryHead2
	character (len=512) :: ModName1, ModName2
	integer :: N1, i

!	write(*,*)'GetReadHeader. NameToRead=', trim(NameToRead)

! ViewHAR does not allow "_" or presumably "@" in header names when
!   edit the file interactively. So drop these.
	ModName1 = NameToRead
	call OmitAllOfSpecifiedString('_', ModName1)
	call OmitAllOfSpecifiedString('@', ModName1)

	if(LEN_TRIM(ModName1).le.llhd) then
	  TryHead1 = ModName1
	else

! Try omitting vowels from the right-hand end until get name of right length
! Do not omit the first letter if it is a vowel
	  do while(LEN_TRIM(ModName1).gt.4.and.LastVowel(ModName1).gt.1)
	    n1=LastVowel(ModName1)
	    if(n1.gt.0) then
	      if(n1.gt.1) then
	        ModName1=ModName1(1:n1-1)//ModName1(n1+1:)
	      end if
	      ModName1 = trim(ModName1)
	    end if
!	    write(*,*)'GetReadHeader. ModName1=', trim(ModName1)
	  end do
	  TryHead1 = ModName1
	end if

! Use TryHead1 if not already used
	if(.not.IsInReadHeaderGPList(TryHead1)) then
	  NewHeader = TryHead1
! Add this to ReadHeaderGPList
	  call AddToReadHeaderGPList( NewHeader)
	  Return
	end if

! Here have not found a new one
! Replace characters at end of TryHead1 until get a new header
! Omit last character from TryHead1 and add digits at the end
	TryHead1 = TryHead1(1:3)
	do i=1,9
	  TryHead1 = trim(TryHead1)//trim(PosinttoString(i))
! Use TryHead1 if not already used
	  if(.not.IsInReadHeaderGPList(TryHead1)) then
	    NewHeader = TryHead1
! Add this to ReadHeaderGPList
	    call AddToReadHeaderGPList( NewHeader)
	    Return
	  end if
	end do

! For the present, return with fatal error
	call InternalErrorStop('GetReadHeader. Did not find new header.')

	return
	end subroutine GetReadHeader

	character (len=255) function ConstructGPName( CoeffVarType, UserSystemType, & 
		Args, SuggestedName, StartIfChangeName)
! Input
! CoeffVarType is 'c' coeff, 'v' variable etc
! UserSystemType is 'u' (user name) or 's' (system name)
	character (len=*) CoeffVarType, UserSystemType, SuggestedName, StartIfChangeName
! Args
	type (domain), pointer :: Args      ! Sets over which arguments range
	

! Use type to work out max length. May change from SuggestedName if it is too long
! Will also change from suggested name if that name has already been used.
! Local
	character (len=LengthGPDeclarationName) TrialName, TrialName2
	type (GPName), pointer :: ThisGPName
	type (domain), pointer :: FoundArgs
	logical :: Found, Finished
	character :: FoundUserSystemType
	integer n1,n2,n3,n4,n5, i
!
! With system names, if SuggestedName is too long for the Type, make up another TrialName
! With user names, must not change them. Is error if name is too long
!
!	write(*,*) 'ConstructGPName: SuggestedName = ', SuggestedName
	Finished = .false.

	if(UserSystemType.eq.'u') then
! Do not change name here, but truncate name if necessary
	  TrialName=SuggestedName
	  ConstructGPName = TrialName
	  Finished = .true.
	  goto 120
	end if
	
! 1. Try omitting ?()"" from TrialName. Then try it
	TrialName = SuggestedName
10	continue
	n1=INDEX(TrialName, '(')
	n2=INDEX(TrialName, ')')
	n3=INDEX(TrialName, '"')
	n4=INDEX(TrialName, '"')
	n5=INDEX(TrialName, '?')
	n1=MAX(n1,n2,n3,n4,n5)
	if(n1.gt.1) then
	  TrialName = TrialName(1:n1-1)//TrialName(n1+1:)
! Try this again
	  goto 10
	end if
! See if TrialName is already one of the names in the list
	call IsNameInGPNameList(TrialName, Found, FoundUserSystemType, FoundArgs)
	if(.not.Found) then
	  ConstructGPName = TrialName
	  Finished = .true.
	end if

	if(.not.Finished) then
! 2. Try replacing ("text") by ?, then replacing ? by integers until get a new one
	  TrialName = SuggestedName
20	  continue
	  n1=INDEX(TrialName, '("')
	  if(n1.gt.1) then
	    n2=INDEX(TrialName(n1+1:),'")')
	    if(n2.gt.0) then
! N2 is position in TrialName(n1+1:), which omits first N1 characters of TrialName. 
!  Hence that is position N1+N2 in TrialName 
	      TrialName=TrialName(1:n1-1)//'?'//TrialName(n1+n2+2:)
! Try this again
	      goto 20
	    end if
	  end if
! Are there any "?". If so, expect just one "?" (but omit later ones if present)
25	  continue
	  n1=INDEX(TrialName,'?')
	  n2=INDEX(TrialName,'?',back=.true.)
	  if(n2.gt.n1) then
! Are two "?". Omit the latter one and try again
	    TrialName=TrialName(1:n2-1)//TrialName(n2+1:)
	    goto 25
	  end if
	  if(n1.gt.1) then
	    if(n2.eq.n1) then
	      TrialName2 = TrialName
	      do i=1,99
	       if(.not.Finished) then
	        TrialName=TrialName2(1:n1-1)//trim(posinttostring(i))//TrialName2(n1+1:)
! See if TrialName is already one of the names in the list
	        call IsNameInGPNameList(TrialName, Found, FoundUserSystemType, FoundArgs)
	        if(.not.Found) then
	          ConstructGPName = TrialName
	          Finished = .true.
	        end if
	       end if ! if(.not.Finished) then
	      end do
	    end if
	  end if
	end if ! if(.not.Finished) then
!
	if(.not.Finished) then
	    call ErrorStop( &
	      'ConstructGPName. Constructing a system name.'// &
	      '\n SuggestedName='//trim(SuggestedName)//'. Need to make up a new name')
	end if
!
!   Here have found new name.
!
	write(*,*)'ConstructGPName. Name used is: ',trim(ConstructGPName)

! Jump here in the user case
120	continue

!  Add this name to those already in GPNameList if it is not already there
	call IsNameInGPNameList(ConstructGPName, Found, FoundUserSystemType, FoundArgs)
	if(.not.Found) then
	  call AddToGPNameList( ConstructGPName, UserSystemType, CoeffVarType, Args)
	end if
!
	return
	end function ConstructGPName

	subroutine IsNameInGPNameList(TrialName, Found, FoundUserSystemType, FoundArgs)
! Input
	character (len=*) TrialName
! Output
	logical :: Found
	character (len=*) :: FoundUserSystemType  ! Only set if Found is true
	type (domain), pointer :: FoundArgs
! Local
	type (GPName), pointer :: ThisGPName
	character (len=LengthGPDeclarationName) TrialName2, ThisName
!
	Found = .false.
	FoundUserSystemType = ' '
	nullify(FoundArgs)

	ThisGPName => GPNameList

! trim TrialName and also shorten it to max allowed and put in lower case for
!		string comparison
	TrialName2 = trim(TrialName)

	do while(associated(ThisGPName).and. .not.Found)
	  ThisName = ThisGPName%name
	  ThisName = trim(ThisName)
	  if(EqualStringsIgnoreCase(TrialName2, ThisName)) then
	    Found = .true.
	    FoundUserSystemType = ThisGPName%UserSystemType
	    FoundArgs => ThisGPName%arg
	  end if
	  ThisGPName => ThisGPName%next
	end do

	end subroutine IsNameInGPNameList

	subroutine AddToGPReadList( Name)
! This adds at the end of GPReadList a GPRead with the specified
!		Name
! Input
	character (len=*) Name
!??	type (domain), pointer :: Args      ! Sets over which arguments range
! Local
	type (GPRead), pointer :: ThisGPRead

	allocate(ThisGPRead)
	ThisGPRead%name = Name
	nullify(ThisGPRead%next)
! Complete info in ThisGPRead
!??	ThisGPRead%arg => Args

	if(.not.associated(GPReadList)) then
	  GPReadList => ThisGPRead
	  GPReadLast => ThisGPRead
	else
	  GPReadLast%next => ThisGPRead
	  GPReadLast => ThisGPRead
	end if
	nullify(GPReadLast%next)

	return
	end subroutine AddToGPReadList

	subroutine IsNameInGPReadList(TrialName, Found)
! Input
	character (len=*) TrialName
! Output
	logical :: Found
! Local
	type (GPRead), pointer :: ThisGPRead
	character (len=LengthGPDeclarationName) TrialName2, ThisName
!
	Found = .false.

	ThisGPRead => GPReadList

! trim TrialName and also shorten it to max allowed and put in lower case for
!		string comparison
	TrialName2 = trim(TrialName)

	do while(associated(ThisGPRead).and. .not.Found)
	  ThisName = ThisGPRead%name
	  ThisName = trim(ThisName)
	  call lower(ThisName)
	  if(ThisName.eq.TrialName2) then
	    Found = .true.
	  end if
	  ThisGPRead => ThisGPRead%next
	end do

	end subroutine IsNameInGPReadList

	subroutine ConvertElastToGPCoeffReal(SOrTType, FunctionString, &
			ElastField, ThisElasticity, ThisElasticityIsConstant)
! Input
	character (len=*) :: SOrTType  ! 's' or 't'  (only used in labelling info)
	character (len=*) :: FunctionString   ! eg "$prod:y(s)" (only used in labelling info)
	type (field), pointer :: ElastField
! Output
	type (GPCoeffVarReal) :: ThisElasticity
	logical :: ThisElasticityIsConstant
! Local
	character (len=255) ElastText
	type (symbol), pointer :: ElastSymbol
!
	ElastText = trim(ElastField%data)

	if(IsReal(ElastText)) then

	  call InitializeGPCoeffVarReal('r', ThisElasticity)
	  ThisElasticity%name = ElastText
	  ThisElasticityIsConstant = .true.

!	 End "if(IsReal(ElastText)) then"
	else

	  call InitializeGPCoeffVarReal('c', ThisElasticity)
	  if(associated(ElastField%ident)) then
	    ThisElasticity%name = ElastField%ident%name
	  else
! Next should be unnecessary ??
!	    write(*,*)'ThisElasticity: ElastField%ident not associated'
	    ThisElasticity%name = ElastField%identtext
! Try to add %ident for this field ElastField
	    call findid(trim(ElastField%identtext), ElastSymbol)
	    if(associated(ElastSymbol)) then
		write(*,*)' Elasticity ident added for ', trim(ElastField%identtext),' - KP'
		ElastField%ident => ElastSymbol
	    end if 
	  end if
	  ThisElasticity%arg => ElastField%arg
	  ThisElasticity%labelinfo = SOrTType//' elasticity in '//trim(FunctionString)
	  ThisElasticity%IsParam = .true.
	  ThisElasticityIsConstant = .false.

!	 End "if(IsReal(ElastText)) then..else.."
	end if

	return
	end subroutine ConvertElastToGPCoeffReal

	subroutine AddToFunctionGPList( ThisFunctionGP)
! This adds ThisFunctionGP at the end of FunctionGPList
! Input
	type (FunctionGP), pointer :: ThisFunctionGP

	if(associated(FunctionGPList)) then
	      FunctionGPLast%next => ThisFunctionGP
	else
	      FunctionGPList => ThisFunctionGP
	end if
	FunctionGPLast => ThisFunctionGP 
	nullify(FunctionGPLast%next)

	return
	end subroutine AddToFunctionGPList

	subroutine AddToNetputGPList( ThisNetputGP)
! This adds ThisNetputGP at the end of NetputGPList
! Input
	type (NetputGP), pointer :: ThisnetputGP

	if(associated(NetputGPList)) then
	      netputGPLast%next => ThisnetputGP
	else
	      NetputGPList => ThisNetputGP
	end if
	NetputGPLast => ThisNetputGP 
	nullify(NetputGPLast%next)

	return
	end subroutine AddToNetputGPList

	subroutine AddToTaxGPList( ThisTaxGP)
! This adds ThisTaxGP at the end of TaxGPList
! Input
	type (TaxGP), pointer :: ThisTaxGP

	if(associated(TaxGPList)) then
	      TaxGPLast%next => ThisTaxGP
	else
	      TaxGPList => ThisTaxGP
	end if
	TaxGPLast => ThisTaxGP 
	nullify(TaxGPLast%next)

	return
	end subroutine AddToTaxGPList

	subroutine AddToEndowmentGPList( ThisEndowmentGP)
! This adds ThisEndowmentGP at the end of EndowmentGPList
! Input
	type (EndowmentGP), pointer :: ThisEndowmentGP

	if(associated(EndowmentGPList)) then
	      EndowmentGPLast%next => ThisEndowmentGP
	else
	      EndowmentGPList => ThisEndowmentGP
	end if
	EndowmentGPLast => ThisEndowmentGP 
	nullify(EndowmentGPLast%next)

	return
	end subroutine AddToEndowmentGPList

	character function ProdSDChar(IsOutput)
!  Supply/Demand character for a o: or i: line in a $prod 
! Input
	logical, intent(in) :: IsOutput

! Either 's' (output) and 'd' (input) OR
!        'o' (output) and 'i' (input)
	if(IsOutput) then
	  ProdSDChar = 's'
	else
	  ProdSDChar = 'd'
	end if

	return
	end function ProdSDChar

	subroutine WriteSymbolList
! Local
	type (symbol), pointer :: CurrentSymbol

!	write(*,*)' Calling WriteSymbolLlist'
	CurrentSymbol => idlist
	do while(associated(CurrentSymbol))
!	  write(*,*)'Symbol name:', trim(CurrentSymbol%name), ' Symbol type:', trim(CurrentSymbol%type)
	  CurrentSymbol => CurrentSymbol%next
	end do
	
	return
	end subroutine WriteSymbolList

	subroutine GetQFieldValueCoeffEtc( &
! Input
	  CurrentQField, &
! CurrentLineField is the field from i:, o: or e: at start of line
! It is CurrentIOField in IO case
	  CurrentLineField, & 
! CurrentLineDomain is CurrentNetput%d in IO case
	  CurrentLineDomain, &
! FunctionAbbreviation is variable name after $prod: or $demand: [eg is "y" in "$prod:y(s)"]
! FunctionString is $prod: or $demand: followed by variable plus arguments [eg is "$prod:y(s)"]
	  FunctionAbbreviation, FunctionString, &
! SupplyDemandEndowmentChar is SupplyDemandChar in IO case
! LineChar is InputOutputChar in IO case
	  SupplyDemandEndowmentChar, LineChar, &
! LineHasRField can only be true if are doing e: line
	  LineHasRField, & 
! Output
! QFieldQuantityCoeff is NewNetputGP%IOQuantityCoeff in IO case
! QFieldQuantityVar is NewNetputGP%IOQFieldQuantityVar in IO case
	  QFieldQuantityCoeff, & 
	  QFieldQuantityVar, & 
! QFieldValueCoeff is NewNetputGP%IOQFieldValueCoeff in IO case
! QFieldValueVar is NewNetputGP%IOQFieldValueVar in IO case
	  QFieldValueCoeff, QFieldValueVar, &
! QFieldArgGP is NewNetputGP%IOArgGP in IO case (equals QFieldValeCoeff%arg at present)
	  QFieldArgGP, &      
! QFieldQuantityCoeffOrRealMPSGE is NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE in IO case
	  QFieldQuantityCoeffOrRealMPSGE, &
! QFieldIsExpression tells if QField is an expression
! QFieldExpressionString - the expression in the q: field (if is such)
	  QFieldIsExpression, QFieldExpressionString, &
	  LinePriceAbbrev)
!

! Input
	type (field), pointer :: CurrentQField, CurrentLineField
	type (domain), pointer :: CurrentLineDomain
	character (len=*), intent(in) :: FunctionAbbreviation, &
	  FunctionString, SupplyDemandEndowmentChar, LineChar
! LineHasRField can only be true if are doing e: line
	logical, intent(in) :: LineHasRField
! Output
	type (GPCoeffVarReal), intent(out) :: QFieldQuantityCoeff, QFieldValueCoeff, QFieldQuantityVar 
	type (GPCoeffVarReal), intent(out) :: QFieldQuantityCoeffOrRealMPSGE, QFieldValueVar
	type (domain), pointer :: QFieldArgGP 
	logical, intent(out) :: QFieldIsExpression
	character (len=*), intent(out) :: QFieldExpressionString 
	character (len=*), intent(out) :: LinePriceAbbrev 
	
! Local
	logical :: NeedGEMPACKCoeff, DomainIsFull, LinePriceIsFullDomain
	type (domain), pointer :: LineFieldArg, LinePriceFullArg
	integer :: n1
	character :: UQChar, PCChar
	logical :: IsELine

! This must work for i: or o: line in $prod or for e: line in $demand
	
	write(*,*)'Calling GetQFieldValueCoeffEtc'

! Initialize
	QFieldIsExpression = .false.
	QFieldExpressionString = ''
	IsELine = (CurrentLineField%label.eq.'e').or.(CurrentLineField%label.eq.'E')

	if(CurrentQField%data(1:1).eq.'(') then
! %data should start with '(' and end with ')'. Strip them off.
	  n1=LEN_TRIM(CurrentQField%data)
	  QFieldExpressionString = CurrentQField%data(2:n1-1)
	  QFieldIsExpression = .true.
	end if

	if(CurrentQField%data(1:1).eq.'-') then
! Sometimes have "-coeff" in q: field in e: line
! All of %data is the expression in this case
	  QFieldExpressionString = CurrentQField%data
	  QFieldIsExpression = .true.
	end if

	if(INDEX(CurrentQField%data,'"').gt.0) then
! Contains a quoted element. Treat this as an expression.
	  QFieldExpressionString = CurrentQField%data
	  write(*,*) 'QFieldExpressionString is: ', trim(QFieldExpressionString)
	  QFieldIsExpression = .true.
	end if

! Look at CurrentLineField and work out if the arguments of the price there
!  are ranging over the full sets as in the declaration of this price in the
!   $commodities part of the file
	LineFieldArg => CurrentLineField%arg
	LinePriceFullArg => CurrentLineField%ident%arg
	if(.not.associated(LinePriceFullArg)) then
	  LinePriceIsFullDomain = .true.
	else
	  LinePriceIsFullDomain = subset(LineFieldArg, LinePriceFullArg) .and. &
				subset(LinePriceFullArg, LineFieldArg)
	end if
! Use LinePriceAbbrev in making up names
	if(LinePriceIsFullDomain) then
	  LinePriceAbbrev = CurrentLineField%identtext   ! eg pf
	else
	  LinePriceAbbrev = CurrentLineField%data        ! eg pf(f1) or pf("lab") or pi("lab",s)
	end if

!	write(*,*)'GetQFieldValueCoeffEtc. CurrentLineField is '//trim(CurrentLineField%label)//':'// &
!	  trim(CurrentLineField%data)
!	write(*,*)'LinePriceAbbrev is: ', LinePriceAbbrev
	
! QFieldQuantityCoeff
! [04/04/04] Now associate QFieldQuantityCoeff with the expression in the q: field since that
!  represents a quantity, not a value (at least in the e: line case)
!QUANT indicates changes made in this connection

! Use what the user gives us if that is of full rank. To check this, compare
!		  CurrentQField%arg and CurrentNetput%d
! Otherwise make up a name

!    NeedGEMPACKCoeff tells if need a GEMPACK Coeff which is different from
!		the MPSGE symbol. This is true if the MPSGE symbol is
!		 numeric (eg "2.1"), in which case it is not in the symbol
!		  table idlist, or if the two domains are not the same.

		NeedGEMPACKCoeff = .false.

! If e: line has r: field, need a gempack coeff
!NO		if(IsEline.and.LineHasRField) then
!NO		  NeedGEMPACKCoeff = .true.
		
! If Q field is an expression, need a gempack coeff
		if(QFieldIsExpression) then
		  NeedGEMPACKCoeff = .true.

		else		

! Note by KP. a real constant (eg "2.1" as in "q:2.1") is not added to the symbol table.
!  See the use of the logical function "numeric" in Tom's code.
!		  write(*,*)'Calling numeric, pos 1'
		  if(numeric(CurrentQField%data)) NeedGEMPACKCoeff = .true.
!		  write(*,*)'Called numeric, pos 1'

		end if

		if(.not.NeedGEMPACKCoeff) then

!  Test if the two domains are the same
		  DomainIsFull = .false.
! If no sets in CurrentNetput%d then DomainIsFull must be true
! IO case		  if(.not.associated(CurrentNetput%d)) then
		  if(.not.associated(CurrentLineDomain)) then
		    DomainIsFull = .true.
		  else
! IO case	    if(.not.IsSubdomainIgnoreOrder(CurrentQField%arg, CurrentNetput%d)) then
		    if(.not.IsSubdomainIgnoreOrder(CurrentQField%arg, CurrentLineDomain)) then
!		      This is an error
		      call ErrorStop('Q field domain is not a subset of Line domain')
		    end if
! IO case	    DomainIsFull = IsSubdomainIgnoreOrder(CurrentNetput%d, CurrentQField%arg)
		    DomainIsFull = IsSubdomainIgnoreOrder(CurrentLineDomain, CurrentQField%arg)
		  end if
		  if(.not.DomainIsFull) NeedGEMPACKCoeff = .true.

!		 End "if(.not.NeedGEMPACKCoeff) then"
		end if

! IO case       call InitializeGPCoeffVarReal('c', NewNetputGP%IOQFieldValueCoeff)
!QUANT	        call InitializeGPCoeffVarReal('c', QFieldValueCoeff)
	        call InitializeGPCoeffVarReal('c', QFieldQuantityCoeff)

! Set up UQChar
		if(LineHasRField) then
! Unrationed quantity
		  UQChar = 'U'
		else
		  UQChar = 'Q'
		end if

		if (.not.NeedGEMPACKCoeff) then
!
!!		    NewNetputGP%IOQFieldValueCoeff%name = CurrentQField%ident%name 
! Use call to ConstructGPName to that this name is added to GPNameList  
! IO case	    NewNetputGP%IOQFieldValueCoeff%name = ConstructGPName('c', 'u', CurrentNetput%d, &
!QUANT		    QFieldValueCoeff%name = ConstructGPName('c', 'u', CurrentLineDomain, &
		    QFieldQuantityCoeff%name = ConstructGPName('c', 'u', CurrentLineDomain, &
		      CurrentQField%ident%name, ' ')
! IO case	    NewNetputGP%IOQFieldValueCoeff%arg => CurrentQField%arg
!QUANT		    QFieldValueCoeff%arg => CurrentQField%arg
		    QFieldQuantityCoeff%arg => CurrentQField%arg
!
!		 End "if (.not.NeedGEMPACKCoeff) then"
		else
!
! IO case	  NewNetputGP%IOQFieldValueCoeff%name = ConstructGPName('c', 's', CurrentNetput%d, &
! IO case	    'V'//SupplyDemandChar//'_'//trim(FunctionAbbreviation)// &
! IO case	    '_'//trim(CurrentLineField%ident%name), 'V'//SupplyDemandChar)
! [01/04/04] If have multiple lines with different parts of the same price 
!     (eg i:pf("lab") and i:pf("cap") or i:pf(f1) where f1 ranges over a subset of f)
! put this into SuggestedName. Routine ConstructGPName replaces these suitably in trying
!  to get a suitable new name.
! So use LinePriceAbbrev in place of CurrentLineFile%ident%name to give ConstructGPName flexibility 
!    to handle multiple lines with different parts of the same price 

!QUANT		  QFieldValueCoeff%name = ConstructGPName('c', 's', CurrentLineDomain, &
		  QFieldQuantityCoeff%name = ConstructGPName('c', 's', CurrentLineDomain, &
! [03/04/04] Drop the "V" at the start since the associated variable QFieldValueVar may have this name
! WAS		    'V'//SupplyDemandEndowmentChar//'_'//trim(FunctionAbbreviation)// &
! [04/04/04] Add "Q"/"U" at start
!QUANT		    SupplyDemandEndowmentChar//'_'//trim(FunctionAbbreviation)// &
		    UQChar//SupplyDemandEndowmentChar//'_'//trim(FunctionAbbreviation)// &
! WAS - old	    '_'//trim(CurrentLineField%ident%name), 'V'//SupplyDemandEndowmentChar)
! WAS		    '_'//trim(LinePriceAbbrev), 'V'//SupplyDemandEndowmentChar)
!QUANT		    '_'//trim(LinePriceAbbrev), SupplyDemandEndowmentChar)
		    '_'//trim(LinePriceAbbrev), UQChar//SupplyDemandEndowmentChar)

! IO case	  NewNetputGP%IOQFieldValueCoeff%arg => CurrentNetput%d 
!QUANT		  QFieldValueCoeff%arg => CurrentLineDomain 
		  QFieldQuantityCoeff%arg => CurrentLineDomain 

		  call GiveWarning(trim(CurrentQField%data)// &
			' will not be updated - better to add full rank symbol here and to read it.')
!
!		 End "if (.not.NeedGEMPACKCoeff) then..else.."
		end if
! IO case       call upper(NewNetputGP%IOQFieldValueCoeff%name)   ! Convert name to upper case
!QUANT	        call upper(QFieldValueCoeff%name)   ! Convert name to upper case
	        call upper(QFieldQuantityCoeff%name)   ! Convert name to upper case
! IO case	NewNetputGP%IOQFieldValueCoeff%labelinfo = & 
! IO case		    'Value (levels) of quantity for '//InputOutputChar//':'//trim(CurrentIOField%data)// &
!QUANT		QFieldValueCoeff%labelinfo = & 
!QUANT		    'Value (levels) of quantity for '//LineChar//':'//trim(CurrentLineField%data)// &
!QUANT		    ' in '//trim(FunctionString)
		if(LineHasRField) then
		  QFieldQuantityCoeff%labelinfo = & 
		    'Rationed quantity (level - no economic meaning) in q: field of '//LineChar//':'//trim(CurrentLineField%data)// &
		    ' in '//trim(FunctionString)
		else
		  QFieldQuantityCoeff%labelinfo = & 
		    'Quantity (level) in q: field of '//LineChar//':'//trim(CurrentLineField%data)// &
		    ' in '//trim(FunctionString)
!  or should it be - ie, give the o: field or the q: field 
!		    'Value (levels) of quantity for '//'q:'//trim(CurrentQField%data)// &
		end if

!-----------------------------------------------------------------------------------------------
! [04/04/04] Now do QFieldValueCoeff - as for QFieldQuantityCoeff with slight differences
	        call InitializeGPCoeffVarReal('c', QFieldValueCoeff)
! Now always a system name
		QFieldValueCoeff%name = ConstructGPName('c', 's', CurrentLineDomain, &
! [03/04/04] Drop the "V" at the start since the associated variable QFieldValueVar may have this name
! WAS		    'V'//SupplyDemandEndowmentChar//'_'//trim(FunctionAbbreviation)// &
! [04/04/0] Drop first "_" to stop getting long too quickly
!DROP		    SupplyDemandEndowmentChar//'_'//trim(FunctionAbbreviation)// &
		    SupplyDemandEndowmentChar//trim(FunctionAbbreviation)// &
! WAS - old	    '_'//trim(CurrentLineField%ident%name), 'V'//SupplyDemandEndowmentChar)
! WAS		    '_'//trim(LinePriceAbbrev), 'V'//SupplyDemandEndowmentChar)
		    '_'//trim(LinePriceAbbrev), SupplyDemandEndowmentChar)
		QFieldValueCoeff%arg => QFieldQuantityCoeff%arg 
		QFieldValueCoeff%labelinfo = & 
		    'Value (levels) of quantity for '//LineChar//':'//trim(CurrentLineField%data)// &
		    ' in '//trim(FunctionString)
!-----------------------------------------------------------------------------------------------

! Set NewNetputGP%IOArgGP - same as NewNetputGP%IOQFieldValueCoeff%arg. Use IOArgGP elsewhere in case
!  later we need to make IOArgGP different from IOQFieldValueCoeff%arg. See the note where
!   IOArgGP is defined under type (netputGP).
! [Do not use CurrentNetput%d because, although these are the right sets, they may be in a different order. See
!    the examples where IOArgGP is defined under type (netputGP).]

! IO case	NewNetputGP%IOArgGP => NewNetputGP%IOQFieldValueCoeff%arg 
		QFieldArgGP => QFieldQuantityCoeff%arg 

! QFieldQuantityVar.
! Set up PCChar
		if(IsELine) then
		  PCChar = 'c' ! QFieldQuantityVar and QFieldValueVar are change variables
		else
		  PCChar = 'p' ! QFieldQuantityVar and QFieldValueVar are %-change variables
		end if
 
! [06/04/04] The names of the variables are always derived from the corresponding coefficient
!  by adding p_ or c_ at the start.
! QFieldQuantityVar
	        call InitializeGPCoeffVarReal('v', QFieldQuantityVar)
		QFieldQuantityVar%name = ConstructGPName('v', 's', QFieldArgGP, &
		      PCChar//'_'//QFieldQuantityCoeff%name, ' ')
	        call lower(QFieldQuantityVar%name)   ! Convert name to lower case
		QFieldQuantityVar%arg => QFieldArgGP 
		QFieldQuantityVar%VarType = PCChar 
		if(LineHasRField) then
		    QFieldQuantityVar%labelinfo = & 
		      'change in unrationed quantity (no economic meaning) for '//LineChar//':'//&
		      trim(CurrentLineField%data)//' in '//trim(FunctionString)
		else
		    QFieldQuantityVar%labelinfo = & 
		      '%-change in quantity for '//LineChar//':'//&
		      trim(CurrentLineField%data)//' in '//trim(FunctionString)
		end if
		QFieldQuantityVar%OrigLevel = QFieldQuantityCoeff%name

! QFieldQuantityVar
	        call InitializeGPCoeffVarReal('v', QFieldValueVar)
		QFieldValueVar%name = ConstructGPName('v', 's', QFieldArgGP, &
		      PCChar//'_'//QFieldValueCoeff%name, ' ')
	        call lower(QFieldValueVar%name)   ! Convert name to lower case
		QFieldValueVar%arg => QFieldArgGP 
		QFieldValueVar%VarType = PCChar 
		if(PCChar.eq.'c') then
		    QFieldValueVar%labelinfo = & 
		      'change in value for '//LineChar//':'//&
		      trim(CurrentLineField%data)//' in '//trim(FunctionString)
		else if(PCChar.eq.'p') then
		    QFieldValueVar%labelinfo = & 
		      '%-change in value for '//LineChar//':'//&
		      trim(CurrentLineField%data)//' in '//trim(FunctionString)
		end if
	        call lower(QFieldQuantityVar%name)   ! Convert name to lower case
		QFieldValueVar%OrigLevel = QFieldValueCoeff%name

! QFieldQuantityCoeffOrRealMPSGE - what the user gives as the initial quantity
! [31/03/04] If QFieldIsExpression, do not (cannot) set up QFieldQuantityCoeffOrRealMPSGE

		if(.not.QFieldIsExpression) then

!		  write(*,*)'Calling numeric, pos 2'
		 if(numeric(CurrentQField%data)) then
! IO case         call InitializeGPCoeffVarReal('r', NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE)
! IO case	  NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%name = CurrentQField%data
	          call InitializeGPCoeffVarReal('r', QFieldQuantityCoeffOrRealMPSGE)
		  QFieldQuantityCoeffOrRealMPSGE%name = CurrentQField%data

		 else

! IO case         call InitializeGPCoeffVarReal('c', NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE)
	          call InitializeGPCoeffVarReal('c', QFieldQuantityCoeffOrRealMPSGE)
! IO case	  NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%name = &
		  QFieldQuantityCoeffOrRealMPSGE%name = &
		    ConstructGPName('c', 'u', CurrentQField%arg, CurrentQField%ident%name, ' ')
! IO case	  NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%arg => CurrentQField%arg
		  QFieldQuantityCoeffOrRealMPSGE%arg => CurrentQField%arg
! IO case	  NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%labelinfo = &
! IO case	    'Pre-simulation Value of quantity for '//InputOutputChar//':'//trim(CurrentIOField%data)
		  QFieldQuantityCoeffOrRealMPSGE%labelinfo = &
		    'Pre-simulation quantity for '//LineChar//':'//trim(CurrentLineField%data)
! IO case	  NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%IsParam = .true.

! QFieldQuantityCoeffOrReal is not a parameter if IsELine and LineHasRField
		  if(.not.(IsELine.and.LineHasRField)) then
		    QFieldQuantityCoeffOrRealMPSGE%IsParam = .true.
		  end if ! if(.not.(IsELine.and.LineHasRField)) then

		 end if ! if(numeric(CurrentQField%data)) then..else..

!		 write(*,*)'Called numeric, pos 2'
		end if ! if(.not.QFieldIsExpression) then

	write(*,*)'Ended GetQFieldValueCoeffEtc'
	

	return
	end subroutine GetQFieldValueCoeffEtc

	subroutine AddToSetGPList( SetName)
! This adds at the end of SetGPList a SetGP with the specified name
! Input
	character (len=*) SetName
! Local
	type (SetGP), pointer :: ThisSetGP

	allocate(ThisSetGP)
	ThisSetGP%name = SetName
	nullify(ThisSetGP%next)

	if(.not.associated(SetGPList)) then
	  SetGPList => ThisSetGP
	  SetGPLast => ThisSetGP
	else
	  SetGPLast%next => ThisSetGP
	  SetGPLast => ThisSetGP
	end if
	nullify(SetGPLast%next)

	return
	end subroutine AddToSetGPList

	subroutine BuildSetGPList

! Local
	type(symbol), pointer :: CurrentSymbol
	type(SetGP), pointer :: NewSetGP, SetGPLast

	NumSets = 0
	nullify(SetGPList)
	nullify(SetGPLast)
	CurrentSymbol => idlist
	do while(associated(CurrentSymbol))  
	  if(CurrentSymbol%type.eq.'set') then
	    NumSets = NumSets + 1
	    call AddToSetGPList(CurrentSymbol%name)
	  end if
	  CurrentSymbol => CurrentSymbol%next
	end do ! while(associated(CurrentSymbol))

	return  
	end subroutine BuildSetGPList

	integer function GetSetNumber(ThisSetName)
! Input
	character (len=*) :: ThisSetName
! Local
	type(setGP), pointer :: CurrentSetGP
	integer :: CurrentSetNumber

	CurrentSetGP => SetGPList
	CurrentSetNumber = 0
	do while(associated(CurrentSetGP))
	  CurrentSetNumber = CurrentSetNumber + 1
	  if(EqualStringsIgnoreCase(ThisSetName, CurrentSetGP%name)) then
	    GetSetNumber = CurrentSetNumber
	    return
	  end if

	  CurrentSetGP => CurrentSetGP%next
	end do ! while(associated(CurrentSetGP))
! Here did not find the set
	GetSetNumber = 0
!  Re-initialize 'pre' messages again at start of each function
	call InitializeMessages('pre',0)
	call InternalErrorStop('GetSetNumber. Did not find set '//trim(ThisSetName)//' in SetGPList.')

	return
	end function GetSetNumber
	
	character (len=llst) function GetSetName(ThisSetNumber)
! Input

	integer, intent(in) :: ThisSetNumber
! Local
	type(setGP), pointer :: CurrentSetGP
	integer :: CurrentSetNumber

	CurrentSetGP => SetGPList
	CurrentSetNumber = 0
	do while(associated(CurrentSetGP))
	  CurrentSetNumber = CurrentSetNumber + 1
	  if(CurrentSetNumber.eq.ThisSetNumber) then
	    GetSetName = CurrentSetGP%name
	    return
	  end if

	  CurrentSetGP => CurrentSetGP%next
	end do ! while(associated(CurrentSetGP))
! Here did not find the number
	GetSetName = ''
!  Re-initialize 'pre' messages again at start of each function
	call InitializeMessages('pre',0)
	call InternalErrorStop('GetSetName. ThisSetNumber is out of range.')

	return
	end function GetSetName

	character (len=llst) function GetSingletonSetName(QuotedElement)
! Input
	character (len=*) :: QuotedElement
! Local
	integer :: n1
	character (len=llst) :: Name1

	n1=LEN_TRIM(QuotedElement)
! For the present, just add "S" to the element name. Later should check
!   that this is a new set.

	Name1 = 's'//QuotedElement(2:n1-1)
	call upper(Name1)
	GetSingletonSetName = Name1

	return
	end function GetSingletonSetName

	subroutine InferAllSubsets

! Local
	type(functionGP), pointer :: CurrentFunctionGP
	type(netputGP), pointer :: CurrentNetputGP
	type(endowmentGP), pointer :: CurrentEndowmentGP
	type(field), pointer :: CurrentDemand
	integer :: NumNetputGPDone, NumEndowmentGPDone
	type(domain), pointer :: IOFieldPriceArg
	type(symbol), pointer :: IOFieldPriceSymbol

! Allocate IsSubsetGP array
	allocate(IsSubsetGP(NumSets,NumSets))
! Initialise it all to false
	IsSubsetGP = .false.

	CurrentFunctionGP => FunctionGPList
	do while(associated(CurrentFunctionGP))

	  if(CurrentFunctionGP%IsProdFunctionGP) then
! Run through the o: and i: lines
	    CurrentNetputGP => CurrentFunctionGP%FirstNetputGP
	    NumNetputGPDone = 0
!	    write(*,*)'InferAllSubsets. Starting NetputGPs'
	    do while(NumNetputGPDone.lt.CurrentFunctionGP%NumIFields+CurrentFunctionGP%NumOFields)
	      IOFieldPriceArg => CurrentNetputGP%netputMPSGE%f%arg
	      IOFieldPriceSymbol => CurrentNetputGP%netputMPSGE%f%ident
! Check arguments against those used when IOFieldPriceName was declared
	      call InferSubsetsGP(IOFieldPriceArg,IOFieldPriceSymbol)
	
	      NumNetputGPDone = NumNetputGPDone + 1
	      CurrentNetputGP => CurrentNetputGP%next
	    end do ! while(NumNetputGPDone.lt.CurrentFunctionGP%NumIFields+CurrentFunctionGP%NumOFields)
	  else if(CurrentFunctionGP%IsDemandFunctionGP) then

	  end if

	  CurrentFunctionGP => CurrentFunctionGP%next
	end do !while(associated(CurrentFunction))

! Run through the functions, working out any implied subsets
	
	end subroutine InferAllSubsets

	subroutine InferSubsetsGP(TheseArgs, ThisSymbol)
! Input
	type(domain), pointer :: TheseArgs
	type(symbol), pointer :: ThisSymbol

! May be able to infer that Sets in TheseArgs must be a subset of the corresponding 
!  set in the arguments of ThisSymbol.
! For example, 
! $sectors
! p(o)
! $prod:y(s)
! o:p(o1) <etc>
! Here can infer that o1 must be a subset of o.

! Local
	integer :: n1,n2,i,ArgNumber, SmallSetNumber, BigSetNumber
	character (len=llst) :: SmallSetName, BigSetName
	type(domain), pointer :: SmallDomain, BigDomain

! Should check that number of sets in TheseArgs is equal to dimension of ThisSymbol
	n1 = GetArgDim(TheseArgs)
	n2 = GetArgDim(ThisSymbol%arg)
	if(n1.ne.n2) then
!  Re-initialize 'pre' messages again at start of each function
	  call InitializeMessages('pre',0)
	  call InternalErrorStop('Mismatch in number of arguments of '//trim(ThisSymbol%name))
	end if
!	write(*,*)'InferSubsetsGP. n1=', n1 

	SmallDomain => TheseArgs
	BigDomain => ThisSymbol%arg
	ArgNumber = 1
	do while(ArgNumber.le.n1)
	  SmallSetName = SmallDomain%s%name
	  BigSetName = BigDomain%s%name
!	  write(*,*)'InferSubsetsGP. SmallSetName=',trim(SmallSetName),' BigSetName=', trim(BigSetName) 
	  if(.not.EqualStringsIgnoreCase(SmallSetName, BigSetName)) then
!		SmallSetName must be a subset of BigSetName
	    SmallSetNumber = GetSetNumber(SmallSetName)
	    BigSetNumber = GetSetNumber(BigSetName)
!	    write(*,*)'InferSubsetsGP. SmallSetName=',trim(SmallSetName),'is subset of BigSetName=', trim(BigSetName) 
	    IsSubsetGP(SmallSetNumber,BigSetNumber) = .true.
	  end if

	  ArgNumber = ArgNumber + 1
	  SmallDomain => SmallDomain%next	  
	  BigDomain => BigDomain%next	  
	end do ! while(ArgNumber.le.n1)

	return
	end subroutine InferSubsetsGP

	subroutine SetupGEMPACKNames
! Sets up names for the various GEMPACK Coefficients and Variables
! Runs through all $prod and all $demand blocks
! Does not make any declarations, just sets up names
	type (field), pointer :: fnf, SElastField, TElastField
	logical :: SElastPresent, TElastPresent, OtherFieldsPresent
	logical :: SElastIsConstant, TElastIsConstant
	type (record), pointer :: r
	type (function), pointer :: fnMPSGE
	type (functionGP), pointer :: NewFnGP
	type (netputGP), pointer :: NewNetputGP
	type (taxGP), pointer :: NewTaxGP
	type (nest), pointer :: n
	type (netputlist), pointer :: e
	type (nest), pointer :: inputs, outputs
	type (domain), pointer :: InputsSubDomain
	type (domain), pointer :: TempDomain

! for $prod
	type (GPCoeffVarReal) :: InputPrice, InputQuantity 
	type (GPCoeffVarReal) :: OutputPrice, OutputQuantity 
	type (GPCoeffVarReal) :: ThisProdVar, ThisTopVar
	type (GPCoeffVarReal) :: ThisElasticity
	type (symbol), pointer :: ElasticitySymbol
! for $demand
	type (GPCoeffVarReal) :: EndowQFieldCoeff, EndowQFieldVar  ! relate to level of q: field
	character :: SupplyDemandEndowmentChar
	character (len=128) LinePriceAbbrev 

	type (netputlist), pointer :: CurrentElements
	type (netput), pointer :: CurrentNetPut
	type (field), pointer :: CurrentIOField, CurrentQField, CurrentEndowField
	type (domain), pointer :: CurrentIODomain, CurrentIOSubdomain
	type (taxlist), pointer :: CurrentTaxList, FirstTaxList, CurrentTaxList2
	logical :: AreSomeTaxes, ThisElasticityIsConstant, DomainIsFull
	logical :: NeedGEMPACKCoeff
	character :: InputOutputChar  ! 'i' for inputs, 'o' for outputs
	character :: SupplyDemandChar  ! 'd' for inputs, 's' for outputs

	type (endowmentlist), pointer :: CurrentEndowmentList
	type (endowment), pointer :: CurrentEndowment
	type (domain), pointer :: CurrentEndowmentDomain
	type (endowmentGP), pointer :: NewEndowmentGP

	character (len=255) :: FunctionAbbreviation  ! [eg is "y" in "$prod:y(s)"]
	type (domain), pointer :: FunctionDomain     ! [eg is "s" in "$prod:y(s)"]
	character (len=255) :: FunctionString  ! [eg is "$prod:y(s)"]
	logical :: IsProd, IsDemand, IsDField, IsEField

	character (len=255) :: CurrentTaxAgentName, CurrentTaxAgentNameGP, ExtraString
	character (len=255) :: TaxAbbrev, TaxLabelPart

	integer :: kpass, i, NumDuplicates

	integer :: fnindex = 0
	character (len=1024) :: rec
	logical DoThis

10015   format(1x,a)     ! KP
10017   format(2x,3a)

	write(*,'(/,1x,a)') 'Calling SetupGEMPACKNames'

	call InitializeMessages('all',0)
	Messages%PostLine(1) = 'This occurred while translating to GEMPACK.'
	Messages%PostLine(2) = 'See the accompanying documentation for more information.'

! Find how many sets there are and build SetGPList
	call BuildSetGPList
!	write(*,*)'Ended BuildSetGPList'

	call WriteSymbolList

	nullify(FunctionGPList)
	nullify(NetputGPList)
	nullify(TaxGPList)
	nullify(EndowmentGPList)
	nullify(GPNameList)
	nullify(GPReadList)
	nullify(GPDeclarationList)
	nullify(ReadHeaderGPList)

	fnMPSGE => functionlist
	do while (associated(fnMPSGE))
	  r => fnMPSGE%firstrecord
	  fnf => r%firstfield

	  IsProd = fnf%label.eq.'$prod'
	  IsDemand = fnf%label.eq.'$demand'

!  Re-initialize 'pre' messages again at start of each function
	  call InitializeMessages('pre',0)
	  Messages%PreLine(1) = 'Processing function '//trim(fnf%label)//':'//trim(fnf%data)

!	  write(*,*)'fnf%label is: ', fnf%label

	  FunctionString = trim(fnf%label)//':'//trim(fnf%data)

	  if (IsProd.or.IsDemand) then

!	    write(*,10015) ' ** SetupGEMPACKNames: Next $prod or $demand'
!	    write(*,10017) 'SetupGEMPACKNames: %data: ', trim(fnf%data)
!	    write(*,10017) 'SetupGEMPACKNames: %identtext: ', trim(fnf%identtext)
!	    write(*,10017) 'SetupGEMPACKNames: %argtext: ', trim(fnf%argtext)


! Get the elasticity info
	    call GetElasticityInfo(fnMPSGE, SElastPresent, TElastPresent, &
		SElastField, TElastField, OtherFieldsPresent)

! OtherFieldsPresent tells if are any fields besides s: and t: fields
	    if(OtherFieldsPresent) then

! This means that have nested function. This is not allowed when translating to GEMPACK.
!  When translating to GEMPACK, first need to run a preliminary program to split the
!   nests into separate functions.

	      call ErrorStop( &
		'Nests not allowed here when translating to GEMPACK.'// &
! "\n" means new line
	        '\n First convert MPS file to split nests into separate functions, then process here again.')

	    end if	    

!	    if(SElastPresent) then
!	      write(*,*)'SElast is present'
!	    else
!	      write(*,*)'SElast is not present'
!	    end if
!	    if(TElastPresent) then
!	      write(*,*)'TElast is present'
!	    else
!	      write(*,*)'TElast is not present'
!	    end if

! Put the info about this into a new FunctionGP
	    call newfunctionGP(NewFnGP, fnMPSGE)

	    NewFnGP%IsProdFunctionGP = IsProd
	    NewFnGP%IsDemandFunctionGP = IsDemand

	    NewFnGP%SElastPresent = SElastPresent
	    NewFnGP%TElastPresent = TElastPresent
	    if(SElastPresent) then
	      call ConvertElastToGPCoeffReal('s', FunctionString, &
		    SElastField, ThisElasticity, ThisElasticityIsConstant)
	      NewFnGP%SElasticity = ThisElasticity
	      NewFnGP%SElastIsConstant = ThisElasticityIsConstant
! Add name to GPNameList
	      call AddToGPNameList( ThisElasticity%name, 'u', 'c', ThisElasticity%arg)
	    end if
	    if(TElastPresent) then
	      call ConvertElastToGPCoeffReal('t', FunctionString, &
		    TElastField, ThisElasticity, ThisElasticityIsConstant)
	      NewFnGP%TElasticity = ThisElasticity
	      NewFnGP%TElastIsConstant = ThisElasticityIsConstant
! Add name to GPNameList
	      call AddToGPNameList( ThisElasticity%name, 'u', 'c', ThisElasticity%arg)
	    end if

! Initialise Count of number of i: and o: fields associated with this function
	    NewFnGP%NumIFields = 0
	    NewFnGP%NumOFields = 0
! Initialise Count of number of d: and e: fields associated with this function
	    NewFnGP%NumDFields = 0
	    NewFnGP%NumEFields = 0

!	   End "if (IsProd.or.IsDemand) then"
	  end if

	  if (IsProd) then

! Temporary
	    if(.not.associated(fnf%ident)) then
	      write(*,*)'$prod. fnf%identtext=', trim(fnf%identtext)
	      call InternalErrorStop('$prod:'//trim(FunctionAbbreviation)//' fnf%ident is not associated.')
	    end if

! Example
! $prod:y(s)          s:sigma(s)         
!  o:p(s)            q:supply(s)
!  i:pf(f)           q:factor(f,s)

! ??
! (temp)    call WriteRecord(' ** $prod record:', r)

	    FunctionAbbreviation = fnf%identtext  ! [eg is "y" in "$prod:y(s)"]
	    FunctionAbbreviation = trim(FunctionAbbreviation)
	    FunctionDomain => fnf%arg    ! [eg is "s" in "$prod:y(s)"]


! ActivityLevelVar   [eg "y(s) in "$prod:y(s)]
!		  ActivityLevelVar should be in the sectors list, so don't
!			to declare it here (nor to set up labelinfo) 
	    call InitializeGPCoeffVarReal('v', NewFnGP%ActivityLevelVar)
	    NewFnGP%ActivityLevelVar%name = trim(FunctionAbbreviation)
	    NewFnGP%ActivityLevelVar%arg => fnf%arg
	    NewFnGP%ActivityLevelVar%labelinfo = fnf%ident%labelinfo

! TotalRevCoeff
	    call InitializeGPCoeffVarReal('c', NewFnGP%TotalRevCoeff)
	    NewFnGP%TotalRevCoeff%name = ConstructGPName('c', 's', FunctionDomain, &
					'R_'//trim(FunctionAbbreviation), 'MR')
	    call upper(NewFnGP%TotalRevCoeff%name)   ! Convert name to upper case
	    NewFnGP%TotalRevCoeff%arg => FunctionDomain 
	    NewFnGP%TotalRevCoeff%labelinfo = 'Total revenue for '//trim(FunctionString)
! MargRevVar
	    call InitializeGPCoeffVarReal('v', NewFnGP%MargRevVar)
	    NewFnGP%MargRevVar%name = ConstructGPName('v', 's', FunctionDomain, &
						'mr_'//trim(FunctionAbbreviation), 'mr')
	    call lower(NewFnGP%MargRevVar%name)   ! Convert name to lower case
	    NewFnGP%MargRevVar%arg => FunctionDomain 
	    NewFnGP%MargRevVar%labelinfo = '%-change in marginal revenue for '//trim(FunctionString)
! TotalCostCoeff
	    call InitializeGPCoeffVarReal('c', NewFnGP%TotalCostCoeff)
	    NewFnGP%TotalCostCoeff%name = ConstructGPName('c', 's', FunctionDomain, &
					'C_'//trim(FunctionAbbreviation), 'MC')
	    call upper(NewFnGP%TotalCostCoeff%name)   ! Convert name to upper case
	    NewFnGP%TotalCostCoeff%arg => FunctionDomain 
	    NewFnGP%TotalCostCoeff%labelinfo = 'Total cost for '//trim(FunctionString)
! MargCostVar
	    call InitializeGPCoeffVarReal('v', NewFnGP%MargCostVar)
	    NewFnGP%MargCostVar%name = ConstructGPName('v', 's', FunctionDomain, &
						'mc_'//trim(FunctionAbbreviation), 'mc')
	    call lower(NewFnGP%MargCostVar%name)   ! Convert name to lower case
	    NewFnGP%MargCostVar%arg => FunctionDomain 
	    NewFnGP%MargCostVar%labelinfo = '%-change in marginal cost for '//trim(FunctionString)

!		Add this to FunctionGPList 
	    call AddToFunctionGPList( NewFnGP)

! Now go through its O: and I: lines adding coeff and var ($prod)
	    do kpass = 1,2
! Do outputs on pass 1, inputs on pass 2
	      if(kpass.eq.1) then
	        CurrentElements => fnMPSGE%outputs%elements
		InputOutputChar = 'o'
	      else if(kpass.eq.2) then
	        CurrentElements => fnMPSGE%inputs%elements
		InputOutputChar = 'i'
	      end if

!	      write(*,*) ' ** KPASS=', kpass, 'InputOutputChar=', InputOutputChar

! Run through the netputs of this type
	      do while(associated(CurrentElements))
	        CurrentNetput => CurrentElements%np
		CurrentIOField => CurrentNetput%f
		CurrentQField => CurrentNetput%qref
		CurrentIODomain => CurrentNetput%d
		CurrentIOSubdomain => CurrentNetput%sd
		CurrentTaxList => CurrentNetput%t
		AreSomeTaxes = associated(CurrentTaxList)

	        Messages%PreLine(2) = 'Processing line beginning '// &
			InputOutputChar//':'//trim(CurrentIOField%data)
!  Re-initialize 'pre' messages past number 2
	        call InitializeMessages('pre',2)

!		write(*,*) 'CurrentNetput%qref is ', trim(CurrentNetput%qref%label), &
!							':', trim(CurrentNetput%qref%data)
!		if(AreSomeTaxes) then
!		  write(*,*)'Are some taxes here'
!		else
!		  write(*,*)'Are no taxes here'
!		end if

!                   SupplyDemandChar is 's' for output, 'd' for input
	        SupplyDemandChar = ProdSDChar(CurrentNetput%output)

! Increment NumIFields or NumOFields
		if(InputOutputChar.eq.'o') NewFnGP%NumOFields = NewFnGP%NumOFields + 1
		if(InputOutputChar.eq.'i') NewFnGP%NumIFields = NewFnGP%NumIFields + 1

!  Begin a new NetputGP
		allocate(NewNetputGP)
		NewNetputGP%netputMPSGE => CurrentNetput

! i: or o: field
! IOPriceVar should be in the commodities list, so don't
!			to declare it here (nor to set up labelinfo) 
	        call InitializeGPCoeffVarReal('v', NewNetputGP%IOPriceVar)
! NO. [This has already been declared as a $commodity] Use call to ConstructGPName to that this name is added to GPNameList  
! NO.		NewNetputGP%IOPriceVar%name = ConstructGPName('v', 'u', CurrentIOField%arg, &
! NO.		    CurrentIOField%ident%name, ' ')
		newNetputGP%IOPriceVar%name = CurrentIOField%ident%name 
	        call lower(NewNetputGP%IOPriceVar%name)   ! Convert name to lower case
		NewNetputGP%IOPriceVar%arg => CurrentIOField%arg
		NewNetputGP%IOPriceVar%labelinfo = CurrentIOField%ident%labelinfo

! IOPriceCoeff
! [This is declared in the $commodities section also]
	        call InitializeGPCoeffVarReal('c', NewNetputGP%IOPriceCoeff)
		newNetputGP%IOPriceCoeff%name = trim(CurrentIOField%ident%name)//'_L' 
	        call upper(NewNetputGP%IOPriceCoeff%name)   ! Convert name to upper case
		NewNetputGP%IOPriceCoeff%arg => CurrentIOField%arg
		NewNetputGP%IOPriceCoeff%labelinfo = CurrentIOField%ident%labelinfo

! q: field
	        Messages%PreLine(3) = 'Processing "'// &
			trim(CurrentQField%label)//':'//trim(CurrentQField%data)//'"'
!  Re-initialize 'pre' messages past number 3
	        call InitializeMessages('pre',3)

! IOQFieldValueCoeff, IOQFieldQuantityVar, QFieldQuantityCoeffOrRealMPSGE
! [28/03/04] Get these via call to routine GetQFieldValueCoeffEtc

		call GetQFieldValueCoeffEtc( &
! Input
	  	  CurrentQField, &
! CurrentLineField is the field from i:, o: or e: at start of line
! It is CurrentIOField in IO case
	  	  CurrentIOField, & 
! CurrentLineDomain is CurrentNetput%d in IO case
	  	  CurrentNetput%d, &
! FunctionAbbreviation is 
	  	  FunctionAbbreviation, FunctionString, &
! SupplyDemandEndowmentChar is SupplyDemandChar in IO case
! LineChar is InputOutputChar in IO case
	  	  SupplyDemandChar, InputOutputChar, &
! LineHasRField must be false here
		  .false., & 
! Output
 ! QFieldQuantityCoeff is NewNetputGP%QFieldIOQuantiyCoeff in IO case
	  	  newNetputGP%IOQFieldQuantityCoeff, & 
 ! QFieldQuantityVar is NewNetputGP%IOQuantiyVar in IO case
	  	  NewNetputGP%IOQFieldQuantityVar, & 
! QFieldValueCoeff is NewNetputGP%IOQFieldValueCoeff in IO case
	  	  newNetputGP%IOQFieldValueCoeff, &
! QFieldValueVar is NewNetputGP%IOQFieldValueVar in IO case
	  	  newNetputGP%IOQFieldValueVar, &
! QFieldArgGP is NewNetputGP%IOArgGP in IO case (equals QFieldValeCoeff%arg at present)
	  	  NewNetputGP%IOArgGP, &      
 ! QFieldQuantityCoeffOrRealMPSGE is NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE in IO case
	  	  NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE, &
		  NewNetputGP%QFieldIsExpression, NewNetputGP%QFieldExpressionString, &
		  LinePriceAbbrev)

!		  if(NewNetputGP%QFieldIsExpression) then
!		    write(*,*) Messages%PreLine(3), 'IsExpression is true'
!		  else
!		    write(*,*) Messages%PreLine(3), 'IsExpression is false'
!		  end if 

!  Add this NetputGP to NetputGPList
		call AddToNetputGPList(NewNetputGP)

! If this is the first netputGP associated with this function, set
		if(NewFnGP%NumIFields+NewFnGP%NumOFields.eq.1) then
		  NewFnGP%FirstNetputGP => NewNetputGP
		end if

! Taxes for this netput

		NewNetputGP%NumTaxTypes = 0
! Record FirstTaxList for use in checking duplicate agent names
		FirstTaxList => CurrentNetput%t

		if(AreSomeTaxes.and. .not.associated(CurrentTaxList)) then
		  call InternalErrorStop('AreSomeTaxes problem')
		end if

		do while(associated(CurrentTaxList))
		  NewNetputGP%NumTaxTypes = NewNetputGP%NumTaxTypes + 1

! At present this GEMPACK conversion cannot handle M: or N: fields in tax stuff
		  if(associated(CurrentTaxList%m).or.associated(CurrentTaxList%n)) then
		    call ErrorStop('Sorry, ths program cannot handle m: or n: tax fields at present.')
		  end if
! Details of this tax type
! For GEMPACK use, can ignore the p: field (CurrentNetput%pref)


! Loop through, processing the pairs of a: and t: fields
!  Use NewTaxGP to hold this info for the current pair
		  allocate(NewTaxGP)
		  NewTaxGP%taxlistMPSGE => CurrentTaxList

! Get agent from a: field. Same agent may be in two or more fields on this line. Add
!  an integer (2,3 etc) to name if so.
		  CurrentTaxAgentName = CurrentTaxList%a%identtext
		  if(NewNetputGP%NumTaxTypes.eq.1) then
		    CurrentTaxAgentNameGP = CurrentTaxAgentName
		  else if(NewNetputGP%NumTaxTypes.gt.1) then
! Go through the earlier ones to see how many times this agent has levied other taxes
		    CurrentTaxList2 => FirstTaxList
		    NumDuplicates = 0	
		    do i=1,NewNetputGP%NumTaxTypes-1
		      if(EqualStringsIgnoreCase(CurrentTaxAgentName, CurrentTaxList2%a%identtext)) then
		        NumDuplicates = NumDuplicates + 1
		      end if
		      if(i.gt.1) then
			CurrentTaxList2 => CurrentTaxList2%next
		      end if
		    end do  ! i
		    if(NumDuplicates.eq.0) then
		      CurrentTaxAgentNameGP = CurrentTaxAgentName
		    else
		      ExtraString = PosintToString(NumDuplicates+1)
		      CurrentTaxAgentNameGP = trim(CurrentTaxAgentName)//ExtraString
		    end if
		  end if

! Details for this tax
!	  type (GPCoeffVarReal) :: TaxPowerCoeff  !! power of tax
!	  type (GPCoeffVarReal) :: TaxRateCoeff  !! rate of tax (number like 0.13)
!	  type (GPCoeffVarReal) :: TaxRevenueCoeff  !! tax revenue
!	  type (GPCoeffVarReal) :: TaxPowerVar  !! %-change in power of tax
!	  type (GPCoeffVarReal) :: TaxRateChangeVar  !! change in tax rate (number like 0.05)
!	  type (GPCoeffVarReal) :: TaxRevenueChangeVar  !! change in tax revenue
! Next will be a real (eg t:0.2) or a Coefficient (eg t:t_p(o))
!	  type (GPCoeffVarReal) :: TaxRateCoeffOrRealMPSGE  !! MPSGE value (initial level) of tax rate

! EXAMPLE. Consider line i:pf(f) q: a:gov t:t_pf(f,s) in $prod:va(s)

! TaxAbbrev. eg is "vapfgov" in line i:pf(f) q: a:gov t:t_pf(f,s) in $prod:va(s)
		  TaxAbbrev = trim(FunctionAbbreviation)//trim(CurrentIOField%identtext)//trim(CurrentTaxAgentNameGP)
! TaxLabelPart. eg is "pf for a:gov t:t_pf(f,s) in $prod:va(s)" in line i:pf(f) q: a:gov t:t_pf(f,s) in $prod:va(s)
		  TaxLabelPart = trim(CurrentIOField%data)//' for "a:'//trim(CurrentTaxAgentName)// &
				   ' t:'//trim(CurrentTaxList%t%data)//'" in '//trim(fnf%label)//':'//trim(fnf%data)
!		  write(*,*)'TaxLabelPart is: ', TaxLabelPart

! NewTaxGP%TaxPowerCoeff.  TIPVAPFGOV
	          call InitializeGPCoeffVarReal('c', NewTaxGP%TaxPowerCoeff)
	          NewTaxGP%TaxPowerCoeff%name = ConstructGPName('c', 's', NewNetputGP%IOArgGP, &
			'T'//InputOutputChar//'P'//trim(TaxAbbrev), 'T'//InputOutputChar//'P')
	          call upper(NewTaxGP%TaxPowerCoeff%name)   ! Convert name to upper case
	          NewTaxGP%TaxPowerCoeff%arg => NewNetputGP%IOArgGP 
	    	  NewTaxGP%TaxPowerCoeff%labelinfo = 'Tax power on '//TaxLabelPart

! NewTaxGP%TaxRateCoeff.  TIRVAPFGOV
	          call InitializeGPCoeffVarReal('c', NewTaxGP%TaxRateCoeff)
! Link name to TaxPowerCoeff, just changing 3rd letter "P" to "R"
	          NewTaxGP%TaxRateCoeff%name = ConstructGPName('c', 's', NewNetputGP%IOArgGP, &
			NewTaxGP%TaxPowerCoeff%name(1:2)//'R'//trim(NewTaxGP%TaxPowerCoeff%name(4:)), 'T'//InputOutputChar//'R')
	          call upper(NewTaxGP%TaxRateCoeff%name)   ! Convert name to upper case
	          NewTaxGP%TaxRateCoeff%arg => NewNetputGP%IOArgGP 
	    	  NewTaxGP%TaxRateCoeff%labelinfo = 'Tax rate on '//TaxLabelPart

! NewTaxGP%TaxRevenueCoeff.  TIVVAPFGOV
	          call InitializeGPCoeffVarReal('c', NewTaxGP%TaxRevenueCoeff)
! Link name to TaxPowerCoeff, just changing 3rd letter "P" to "V"
	          NewTaxGP%TaxRevenueCoeff%name = ConstructGPName('c', 's', NewNetputGP%IOArgGP, &
			NewTaxGP%TaxPowerCoeff%name(1:2)//'V'//trim(NewTaxGP%TaxPowerCoeff%name(4:)), 'T'//InputOutputChar//'V')
	          call upper(NewTaxGP%TaxRevenueCoeff%name)   ! Convert name to upper case
	          NewTaxGP%TaxRevenueCoeff%arg => NewNetputGP%IOArgGP 
	    	  NewTaxGP%TaxRevenueCoeff%labelinfo = 'Tax revenue on '//TaxLabelPart

! NewTaxGP%TaxPowerVar.  tivapfgov
	          call InitializeGPCoeffVarReal('v', NewTaxGP%TaxPowerVar)
	          NewTaxGP%TaxPowerVar%name = ConstructGPName('v', 's', NewNetputGP%IOArgGP, &
! Make this name always related to the corresponding coeff. Drop the P which is 3rd letter
		       NewTaxGP%TaxPowerCoeff%name(1:2)//NewTaxGP%TaxPowerCoeff%name(4:), &
		           NewTaxGP%TaxPowerCoeff%name(1:2))
	          call lower(NewTaxGP%TaxPowerVar%name)   ! Convert name to lower case
	          NewTaxGP%TaxPowerVar%arg => NewNetputGP%IOArgGP 
	    	  NewTaxGP%TaxPowerVar%labelinfo = '%-change in power of tax on '//TaxLabelPart
	    	  NewTaxGP%TaxPowerVar%OrigLevel = trim(NewTaxGP%TaxPowerCoeff%name)  

! NewTaxGP%TaxRateChangeVar.  c_tirvapfgov
	          call InitializeGPCoeffVarReal('v', NewTaxGP%TaxRateChangeVar)
	          NewTaxGP%TaxRateChangeVar%name = ConstructGPName('v', 's', NewNetputGP%IOArgGP, &
! Make this name always related to the corresponding coeff. Just add "c_" at the start.
		       "c_"//trim(NewTaxGP%TaxRateCoeff%name), "c_"//trim(NewTaxGP%TaxRateCoeff%name) )
	          call lower(NewTaxGP%TaxRateChangeVar%name)   ! Convert name to lower case
	          NewTaxGP%TaxRateChangeVar%arg => NewNetputGP%IOArgGP 
	    	  NewTaxGP%TaxRateChangeVar%labelinfo = 'Change in rate of tax on '//TaxLabelPart
	    	  NewTaxGP%TaxRateChangeVar%VarType = 'c'
	    	  NewTaxGP%TaxRateChangeVar%OrigLevel = trim(NewTaxGP%TaxRateCoeff%name)  

! NewTaxGP%TaxRevenueChangeVar.  c_tivvapfgov
	          call InitializeGPCoeffVarReal('v', NewTaxGP%TaxRevenueChangeVar)
	          NewTaxGP%TaxRevenueChangeVar%name = ConstructGPName('v', 's', NewNetputGP%IOArgGP, &
! Make this name always related to the corresponding coeff. Just add "c_" at the start.
		       "c_"//trim(NewTaxGP%TaxRevenueCoeff%name), "c_"//trim(NewTaxGP%TaxRevenueCoeff%name) )
	          call lower(NewTaxGP%TaxRevenueChangeVar%name)   ! Convert name to lower case
	          NewTaxGP%TaxRevenueChangeVar%arg => NewNetputGP%IOArgGP 
	    	  NewTaxGP%TaxRevenueChangeVar%labelinfo = 'Change in revenue from tax on '//TaxLabelPart
	    	  NewTaxGP%TaxRevenueChangeVar%VarType = 'c'
	    	  NewTaxGP%TaxRevenueChangeVar%OrigLevel = trim(NewTaxGP%TaxRevenueCoeff%name)  

! TaxRateCoeffOrRealMPSGE - what the user gives as the initial value
		  if(numeric(CurrentTaxList%t%data)) then
	            call InitializeGPCoeffVarReal('r', NewTaxGP%TaxRateCoeffOrRealMPSGE)
		    NewTaxGP%TaxRateCoeffOrRealMPSGE%name = CurrentTaxList%t%data
		  else
	            call InitializeGPCoeffVarReal('c', NewTaxGP%TaxRateCoeffOrRealMPSGE)
		    NewTaxGP%TaxRateCoeffOrRealMPSGE%name = &
		      ConstructGPName('c', 'u', CurrentTaxList%t%arg, CurrentTaxList%t%ident%name, ' ')
		    NewTaxGP%TaxRateCoeffOrRealMPSGE%arg => CurrentTaxList%t%arg

! labelinfo. This can easily get too long, so drop "value of" if necessary
		    if(LEN_TRIM('Pre-sim value of tax rate on '//TaxLabelPart).le.73) then
		      NewTaxGP%TaxRateCoeffOrRealMPSGE%labelinfo = 'Pre-sim value of tax rate on '//TaxLabelPart
		    else
		      NewTaxGP%TaxRateCoeffOrRealMPSGE%labelinfo = 'Pre-sim tax rate on '//TaxLabelPart
		    end if
		    NewTaxGP%TaxRateCoeffOrRealMPSGE%IsParam = .true.
		  end if

		  if(NewNetputGP%NumTaxTypes.ge.2) then
!
! Next only needed if are 2 or more tax types for this netput
!	  type (GPCoeffVarReal) :: TotalTaxPowerCoeff  !! power of tax (all tax types)
!	  type (GPCoeffVarReal) :: TotalTaxPowerVar  !! %-change in power of tax (all tax types)

! TotalTaxPowerCoeff TIP_VA_PF
	            call InitializeGPCoeffVarReal('c', NewNetputGP%TotalTaxPowerCoeff)
	            NewNetputGP%TotalTaxPowerCoeff%name = ConstructGPName('c', 's', NewNetputGP%IOArgGP, &
		      'T'//InputOutputChar//'P_'//trim(FunctionAbbreviation)//'_'//trim(CurrentIOField%ident%name), &
		      'T'//InputOutputChar//'P')
	            call upper(NewNetputGP%TotalTaxPowerCoeff%name)   ! Convert name to upper case
	            NewNetputGP%TotalTaxPowerCoeff%arg => NewNetputGP%IOArgGP 
	    	    NewNetputGP%TotalTaxPowerCoeff%labelinfo = 'Total power of taxes on '// &
		      trim(CurrentIOField%data)//' in '//trim(fnf%label)//':'//trim(fnf%data)

! TotalTaxPowerVar ti_va_pf
	            call InitializeGPCoeffVarReal('v', NewNetputGP%TotalTaxPowerVar)
	            NewNetputGP%TotalTaxPowerVar%name = ConstructGPName('c', 's', NewNetputGP%IOArgGP, &
		      't'//InputOutputChar//'_'//trim(FunctionAbbreviation)//'_'//trim(CurrentIOField%ident%name), &
		      't'//InputOutputChar)
	            call lower(NewNetputGP%TotalTaxPowerVar%name)   ! Convert name to lower case
	            NewNetputGP%TotalTaxPowerVar%arg => NewNetputGP%IOArgGP 
	    	    NewNetputGP%TotalTaxPowerVar%labelinfo = '%-change in total power of taxes on '// &
		      trim(CurrentIOField%data)//' in '//trim(fnf%label)//':'//trim(fnf%data)
	    	    NewNetputGP%TotalTaxPowerVar%OrigLevel = NewNetputGP%TotalTaxPowerCoeff%name

!		   End "if(NewNetputGP%NumTaxTypes.ge.2) then"
		  end if


! Set NewNetputGP%FirstTaxGP for the current NetputGP
		  if(NewNetputGP%NumTaxTypes.eq.1) then
		    NewNetputGP%FirstTaxGP => NewTaxGP
		  end if

! Add this TaxGP to TaxGPList
		  call AddToTaxGPList(NewTaxGP)

		  CurrentTaxList => CurrentTaxList%next 
		end do  ! while(associated(CurrentTaxList))

		if(NewNetputGP%NumTaxTypes.gt.0) then

!		  write(*,*)'NumTaxtypes for '//InputOutputChar//':'//trim(CurrentIOField%data)//' is ',NewNetputGP%NumTaxTypes 
!
! Next only needed if are one or more taxing agents
!	  type (GPCoeffVarReal) :: IOPriceIncTaxesVar    !! %-change in price including taxes
! Next only needed for outputs when are one or more taxing agents. 
!	  type (GPCoeffVarReal) :: IOQFieldValueVar    !! %-change in value of input or output

! Leave these until now so can use TotalTaxPowerCoeff in OrigLevels below.

! IOPriceIncTaxesVar
! Example. Consider  o:p(o) q:supply(s,o) + taxes in $prod:y(s). 
!   Then IOQFieldQuantityVar%name is p_vs_y_p(s,o) and so we want IOPriceIncTaxesVar%name to be p_ps_y_p(s,o).  
	          call InitializeGPCoeffVarReal('v', NewNetputGP%IOPriceIncTaxesVar)
! Base name on that for IOQFieldValueVar - just change "v" in position 3 to "p"
		  NewNetputGP%IOPriceIncTaxesVar%name = ConstructGPName('v', 's', NewNetputGP%IOArgGP, &
		    'p_p'//trim(NewNetputGp%IOQFieldValueVar%name(4:)), 'p_p'//trim(NewNetputGp%IOQFieldQuantityVar%name(4:)))
		  NewNetputGP%IOPriceIncTaxesVar%arg => NewNetputGP%IOArgGP 
		  NewNetputGP%IOPriceIncTaxesVar%labelinfo = & 
		    '%-change in price inc taxes for '//InputOutputChar//':'//&
		    trim(CurrentIOField%data)//' in '//trim(FunctionString)
	          call lower(NewNetputGP%IOPriceIncTaxesVar%name)   ! Convert name to lower case
		  if(NewNetputGP%NumTaxTypes.eq.1) then
		    NewNetputGP%IOPriceIncTaxesVar%OrigLevel = NewNetputGP%FirstTaxGP%TaxPowerCoeff%name
		  else if(NewNetputGP%NumTaxTypes.ge.2) then
		    NewNetputGP%IOPriceIncTaxesVar%OrigLevel = NewNetputGP%TotalTaxPowerCoeff%name
		  end if  


! [06/04/04] IOQFieldValueVar is now set up in call to GetQFieldValueCoeff so does not need to be done here.
! IOQFieldValueVar (includes taxes in o: field, does not include taxes in i: field)
! Example. Consider  o:p(o) q:supply(s,o) + taxes in $prod:y(s). Then IOQFieldValueVar%name is vs_y_p(s,o)  
!EARLIER	            call InitializeGPCoeffVarReal('v', NewNetputGP%IOQFieldValueVar)
! Base name on that for IOQFieldQuantityVar - just change "q" at start to "v"
!EARLIER		    NewNetputGP%IOQFieldValueVar%name = ConstructGPName('v', 's', NewNetputGP%IOArgGP, &
!EARLIER		      'v'//trim(NewNetputGp%IOQFieldQuantityVar%name(2:)), 'v'//trim(NewNetputGp%IOQFieldQuantityVar%name(2:)))
!EARLIER		    NewNetputGP%IOQFieldValueVar%arg => NewNetputGP%IOArgGP 
!EARLIER		    NewNetputGP%IOQFieldValueVar%labelinfo = & 
!EARLIER		      '%-change in value inc taxes for '//InputOutputChar//':'//&
!EARLIER		      trim(CurrentIOField%data)//' in '//trim(FunctionString)
!EARLIER	            call lower(NewNetputGP%IOQFieldValueVar%name)   ! Convert name to lower case
!EARLIER		    NewNetputGP%IOQFieldValueVar%OrigLevel = NewNetputGP%IOQFieldValueCoeff%name  


!		 End "if(NewNetputGP%NumTaxTypes.gt.0) then"
		end if



		CurrentElements => CurrentElements%next
	      end do

	    end do ! kpass loop

!	   End "if (IsProd) then"
	  else if (IsDemand) then

! (temp)    call WriteRecord(' ** $demand record:', r)

	    FunctionAbbreviation = fnf%identtext  ! [eg is "y" in "$prod:y(s)"]
	    FunctionAbbreviation = trim(FunctionAbbreviation)
	    FunctionDomain => fnf%arg    ! [eg is "s" in "$prod:y(s)"]
!	    write(*,*)'$demand. FunctionAbbreviation is: ', trim(FunctionAbbreviation)
	    if(.not.associated(fnf%ident)) then
	      write(*,*)'$demand. fnf%identtext=', trim(fnf%identtext)
	      write(*,*) '$demand:'//trim(FunctionAbbreviation)//' fnf%ident is not associated. Should stop.'
!??	      call InternalErrorStop('$demand:'//trim(FunctionAbbreviation)//' fnf%ident is not associated.')
	    end if

! Example
! $demand:ra
!  d:pc              q:cons
!  e:pf(f)           q:endow(f)

! IncomeVar   [eg "ra" in "$demand:ra"]
!		  IncomeVar should be in the sectors list, so don't
!			to declare it here (nor to set up labelinfo) 
	    call InitializeGPCoeffVarReal('v', NewFnGP%IncomeVar)
	    NewFnGP%IncomeVar%name = trim(FunctionAbbreviation)
	    NewFnGP%IncomeVar%arg => fnf%arg
	    NewFnGP%IncomeVar%labelinfo = fnf%ident%labelinfo

! I think that there must always be exactly one d: field
	    NewFnGP%NumDFields = 1



! TotalIncomeCoeff
!	  type (GPCoeffVarReal) :: TotalIncomeCoeff !! Total income for this consumer
		! Is equal to total value of endowments (summed over all e: lines) plus total tax income
! One for each $demand function
! The value of this is the sum over ALL e: fields.
! The %-change in this is just IncomeVar. So this could be a levels variable.

	    call InitializeGPCoeffVarReal('c', NewFnGP%TotalIncomeCoeff)
	    NewFnGP%TotalIncomeCoeff%name = ConstructGPName('c', 's', FunctionDomain, &
					'VI_'//trim(FunctionAbbreviation), 'VI')
	    call upper(NewFnGP%TotalIncomeCoeff%name)   ! Convert name to upper case
	    NewFnGP%TotalIncomeCoeff%arg => FunctionDomain 
	    NewFnGP%TotalIncomeCoeff%labelinfo = 'Total value of income for consumer in '//trim(FunctionString)
 


! Other things
! for each e: line, have EndowQFieldQuantityVar eg qs_ra_pf

! In any $demand: block, is just one equation and it is a natural levels equation.
! It says TotalIncomeCoeff = SUM value from q: fields, summed across all e: lines
!                                        PLUS the sum of all revenue accruing to the agent (look in $prod: blocks)

! Note that the prices in $commodities are MARKET prices !!

! Will have linear eq linking %-change in values in q: fields with the price in e: field and the
!   EndowQFieldQuantityVar  


!		Add this to FunctionGPList 
	    call AddToFunctionGPList( NewFnGP)

	    CurrentEndowmentList => fnMPSGE%endowments

! Run through the endowments in this list
	    do while(associated(CurrentEndowmentList))
	        CurrentEndowment => CurrentEndowmentList%e
		CurrentQField => CurrentEndowment%q
		CurrentEndowmentDomain => CurrentEndowment%d
	        CurrentEndowField => CurrentEndowment%f

	        NewFnGP%NumEFields = NewFnGP%NumEFields + 1


!  Begin a new EndowmentGP
		allocate(NewEndowmentGP)
		NewEndowmentGP%endowmentMPSGE => CurrentEndowment

	        Messages%PreLine(2) = 'Processing line beginning '// &
			'e:'//trim(CurrentEndowField%data)
!  Re-initialize 'pre' messages past number 2
	        call InitializeMessages('pre',2)

!	type endowmentGP
! Additions to Tom Rutherford's type endowment. These additions needed for GEMPACK
! This is actually a list of endowmentGP's
!	  type(endowment), pointer :: endowmentMPSGE !! Pointer to the MPSGE endowment

! EndowArgGP is analagous to IOArgGP
! At present, EndowArgGP and EndowQFieldValueCoeff%arg are always the same. But use EndowArgGP in the code
!  in case this changes in some special cases which are excluded at present.
!	  type (domain), pointer :: EndowArgGP  

! The next 4 are analagous to the corresponding IO things in a netputGP
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityVar  !! %-change in quantity of each endowment  
!	  type (GPCoeffVarReal) :: EndowPriceVar    !! %-change in price (The price after e: at start of line)
!	  type (GPCoeffVarReal) :: EndowQFieldValueCoeff  !! value (level) of endowment
!	  type (GPCoeffVarReal) :: EndowQFieldValueVar  !! %-change in value (level) of endowment
! Next will be a real (eg q:200) or a Coefficient (eg q:totfac)
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeffOrRealMPSGE  !! MPSGE value (initial level) of endowment

!	  type(endowmentGP), pointer :: next    !! Pointer to next in linked list of all endowmentGPs
!	end type endowmentGP

! EndowPriceVar is price after e: at start of this e: line
!		  EndowPriceVar should be in the commodities list, so don't need
!			to declare it here (nor to set up labelinfo) 
	        call InitializeGPCoeffVarReal('v', NewEndowmentGP%EndowPriceVar)
		NewEndowmentGP%EndowPriceVar%name = CurrentEndowField%ident%name 
	        call lower(NewEndowmentGP%EndowPriceVar%name)   ! Convert name to lower case
		NewEndowmentGP%EndowPriceVar%arg => CurrentEndowField%arg
		NewEndowmentGP%EndowPriceVar%labelinfo = CurrentEndowField%ident%labelinfo

! EndowPriceCoeff - levels version of price after e: at start of this e: line
!		  EndowPriceCoeff is also declared in the $commodities section so don't need
!			to declare it here (nor to set up labelinfo) 
	        call InitializeGPCoeffVarReal('c', NewEndowmentGP%EndowPriceCoeff)
		NewEndowmentGP%EndowPriceCoeff%name = trim(CurrentEndowField%ident%name)//'_L' 
	        call upper(NewEndowmentGP%EndowPriceCoeff%name)   ! Convert name to upper case
		NewEndowmentGP%EndowPriceCoeff%arg => CurrentEndowField%arg
		NewEndowmentGP%EndowPriceCoeff%labelinfo = CurrentEndowField%ident%labelinfo

! q: field
	        Messages%PreLine(3) = 'Processing "'// &
			trim(CurrentQField%label)//':'//trim(CurrentQField%data)//'"'
!  Re-initialize 'pre' messages past number 3
	        call InitializeMessages('pre',3)

! EndowQFieldValueCoeff, EndowQFieldQuantityVar, EndowQFieldQuantityCoeffOrRealMPSGE
! [28/03/04] Get these via call to routine GetQFieldValueCoeffEtc

! Is there an r: field
		NewEndowmentGP%HasRField = (associated(NewEndowmentGP%EndowmentMPSGE%r))

		SupplyDemandEndowmentChar = 's'

		call GetQFieldValueCoeffEtc( &
! Input
	  	  CurrentQField, &
! CurrentLineField is the field from i:, o: or e: at start of line
! It is CurrentIOField in IO case
	  	  CurrentEndowField, & 
! CurrentLineDomain is CurrentNetput%d in IO case
	  	  CurrentEndowmentDomain, &
! FunctionAbbreviation is variable name after $prod: or $demand: [eg is "y" in "$prod:y(s)"]
! FunctionString is $prod: or $demand: followed by variable plus arguments [eg is "$prod:y(s)"]
! These are already set up in this case 
	  	  FunctionAbbreviation, FunctionString, &
! SupplyDemandEndowmentChar is SupplyDemandChar in IO case
! LineChar is InputOutputChar in IO case
! call	  	  SupplyDemandEndowmentChar, LineChar, & 
	  	  SupplyDemandEndowmentChar, 'e', & 
! LineHasRField (in the call)
		  NewEndowmentGP%HasRField, &
! Output
! EndowQFieldCoeff is either EndowQFieldQuantityCoeff or EndowUnrationedQuantityCoeff
	  	  EndowQFieldCoeff, &
! EndowQFieldVar is either EndowQFieldQuantityVar or EndowUnrationedQuantityVar
	  	  EndowQFieldVar, & 
! QFieldValueCoeff is NewNetputGP%IOQFieldValueCoeff in IO case
	  	  NewEndowmentGP%EndowQFieldValueCoeff, &
! QFieldValueVar is NewNetputGP%IOQFieldValueVar in IO case
	  	  NewEndowmentGP%EndowQFieldValueVar, &
! QFieldArgGP is NewNetputGP%IOArgGP in IO case (equals QFieldValeCoeff%arg at present)
	  	  NewEndowmentGP%EndowArgGP, &      
! QuantityVar is NewNetputGP%IOQFieldQuantityVar in IO case
! QFieldQuantityCoeffOrRealMPSGE is NewNetputGP%IOQFieldQuantityCoeffOrRealMPSGE in IO case
	  	  NewEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE, &
		  NewEndowmentGP%QFieldIsExpression, NewEndowmentGP%QFieldExpressionString, &
		  LinePriceAbbrev)

!		  if(NewNetputGP%QFieldIsExpression) then
!		    write(*,*) Messages%PreLine(3), 'IsExpression is true'
!		  else
!		    write(*,*) Messages%PreLine(3), 'IsExpression is false'
!		  end if 

		if(NewEndowmentGP%HasRField) then
! Q: field is unrationed quantity
	          NewEndowmentGP%EndowUnrationedQuantityCoeff = EndowQFieldCoeff
	          NewEndowmentGP%EndowUnrationedQuantityVar = EndowQFieldVar
		  write(*,*)'EndowUnrationedQuantityCoeff%name is ', NewEndowmentGP%EndowUnrationedQuantityCoeff%name
		  write(*,*)'EndowUnrationedQuantityVar%name is ', NewEndowmentGP%EndowUnrationedQuantityVar%name

! Need to add EndowQFieldQuantityCoeff and EndowQFieldQuantityVar
	          call InitializeGPCoeffVarReal('c', NewEndowmentGP%EndowQFieldQuantityCoeff)
! Always attempt to use the full system name here
	          NewEndowmentGP%EndowQFieldQuantityCoeff%name = ConstructGPName('c', 's', EndowQFieldCoeff%arg, &
! Add "Q" at start
		    'Q'//SupplyDemandEndowmentChar//'_'//trim(FunctionAbbreviation)// &
		    '_'//trim(LinePriceAbbrev), 'Q'//SupplyDemandEndowmentChar)
	          call upper(NewEndowmentGP%EndowQFieldQuantityCoeff%name)   ! Convert name to upper case
	          NewEndowmentGP%EndowQFieldQuantityCoeff%arg => EndowQFieldCoeff%arg 
	          NewEndowmentGP%EndowQFieldQuantityCoeff%labelinfo = 'Actual quantity in e:'// &
			trim(CurrentEndowField%data)//' in '//trim(FunctionString)

	          call InitializeGPCoeffVarReal('v', NewEndowmentGP%EndowQFieldQuantityVar)
	          NewEndowmentGP%EndowQFieldQuantityVar%name = ConstructGPName('v', 's', EndowQFieldCoeff%arg, &
					'c_'//NewEndowmentGP%EndowQFieldQuantityCoeff%name, ' ')
	          call lower(NewEndowmentGP%EndowQFieldQuantityVar%name)   ! Convert name to lower case
	          NewEndowmentGP%EndowQFieldQuantityVar%arg => EndowQFieldCoeff%arg 
	          NewEndowmentGP%EndowQFieldQuantityVar%VarType = 'c'   ! is a change variable
	          NewEndowmentGP%EndowQFieldQuantityVar%labelinfo = 'Change in actual quantity in e:'// &
			trim(CurrentEndowField%data)//' in '//trim(FunctionString)

		else
! Q: field is the actual quantity
	          NewEndowmentGP%EndowQFieldQuantityCoeff = EndowQFieldCoeff
		  NewEndowmentGP%EndowQFieldQuantityVar = EndowQFieldVar
		end if
   
		if(NewEndowmentGP%HasRField) then
! In this case, have the levels equation EndowQFieldQuantityCoeff = EndowUnrationedQuantityCoeff * EndowRFieldCoeff 
		end if

! [06/04/04] EndowQFieldValueVar is now set up in call to GetQFieldValueCoeffEtc so do not need to redo it here
! EndowQFieldValueVar
! Example. Consider  o:p(o) q:supply(s,o) + taxes in $prod:y(s). Then IOQFieldValueVar%name is vs_y_p(s,o)  
!EARLIER	        call InitializeGPCoeffVarReal('v', NewEndowmentGP%EndowQFieldValueVar)
! Base name on that for EndowQFieldQuantityVar - just change "q" at start to "v"
!EARLIER		NewEndowmentGP%EndowQFieldValueVar%name = ConstructGPName('v', 's', NewEndowmentGP%EndowArgGP, &
!EARLIER		      'v'//trim(NewEndowmentGP%EndowQFieldQuantityVar%name(2:)), 'v'//trim(NewEndowmentGP%EndowQFieldQuantityVar%name(2:)))
!EARLIER		NewEndowmentGP%EndowQFieldValueVar%arg => NewEndowmentGP%EndowArgGP 
!EARLIER		NewEndowmentGP%EndowQFieldValueVar%labelinfo = & 
!EARLIER		      '%-change in value of e:'//&
!EARLIER		      trim(CurrentEndowField%data)//' in '//trim(FunctionString)
!EARLIER	        call lower(NewEndowmentGP%EndowQFieldValueVar%name)   ! Convert name to lower case
!EARLIER		NewEndowmentGP%EndowQFieldValueVar%OrigLevel = NewEndowmentGP%EndowQFieldValueCoeff%name

! EndowRFieldCoeff and EndowRFieldVar
		if(NewEndowmentGP%HasRField) then
	          call InitializeGPCoeffVarReal('c', NewEndowmentGP%EndowRFieldCoeff)
	          NewEndowmentGP%EndowRFieldCoeff%name = ConstructGPName('c', 'u', NewEndowmentGP%EndowmentMPSGE%r%arg, &
					NewEndowmentGP%EndowmentMPSGE%r%identtext, ' ')
	          call upper(NewEndowmentGP%EndowRFieldCoeff%name)   ! Convert name to upper case
	          NewEndowmentGP%EndowRFieldCoeff%arg => NewEndowmentGP%EndowmentMPSGE%r%arg 
	          NewEndowmentGP%EndowRFieldCoeff%labelinfo = 'Rationing quantity in e:'// &
			trim(CurrentEndowField%data)//' in '//trim(FunctionString)

	          call InitializeGPCoeffVarReal('v', NewEndowmentGP%EndowRFieldVar)
	          NewEndowmentGP%EndowRFieldVar%name = ConstructGPName('v', 'u', NewEndowmentGP%EndowmentMPSGE%r%arg, &
					'c_'//NewEndowmentGP%EndowmentMPSGE%r%identtext, ' ')
	          call lower(NewEndowmentGP%EndowRFieldVar%name)   ! Convert name to lower case
	          NewEndowmentGP%EndowRFieldVar%arg => NewEndowmentGP%EndowmentMPSGE%r%arg 
	          NewEndowmentGP%EndowRFieldVar%VarType = 'c'   ! is a change variable
	          NewEndowmentGP%EndowRFieldVar%labelinfo = 'Change in rationing quantity in e:'// &
			trim(CurrentEndowField%data)//' in '//trim(FunctionString)
		end if

!  Add this EndowmentGP to EndowmentGPList
		call AddToEndowmentGPList(NewEndowmentGP)

! If this is the first endowmentGP associated with this function, set
		if(NewFnGP%NumEFields.eq.1) then
		  NewFnGP%FirstEndowmentGP => NewEndowmentGP
		end if


	        CurrentEndowmentList => CurrentEndowmentList%next
	    end do 

!	   End "if (IsProd) then..else if(IsDemand) then.."
	  end if

	  fnMPSGE => fnMPSGE%next
	end do

! Finish with WalrasSlackVar
	call InitializeGPCoeffVarReal('v', WalrasSlackVar)
	nullify(WalrasSlackVar%arg)
	WalrasSlackVar%name = ConstructGPName('v', 's', WalrasSlackVar%arg, &
					'walrasslack', ' ')
	call lower(WalrasSlackVar%name)   ! Convert name to lower case
	WalrasSlackVar%labelinfo = 'Check of omitted market clearing equation. Result should be zero.'

! Write details about FunctionGPList and GPNameList

	call WriteFunctionGPList
	call WriteGPNameList

! Infer any subsets - must wait until FunctionGPList is set up to do this
	call InferAllSubsets
!	write(*,*)'Ended InferAllSubsets'

	
	write(*,*)'Ended SetupGEMPACKNames'

	return
	end subroutine SetupGEMPACKNames

	subroutine WriteGPNameList

	type (GPName), pointer :: ThisGPName

	ThisGPName => GPNameList
	do while(associated(ThisGPName))

!	  write(*,*) 'WriteGPNameList: ThisGPName%name is ', ThisGPName%name

	  ThisGPName => ThisGPName%next
	end do

	end subroutine WriteGPNameList

	subroutine WriteFunctionGPList

! Writes some details about everything in FunctionGPList

	type (functionGP), pointer :: fnGP

!	write(*,*)' Begin WriteFunctionGPList'

	fnGP => FunctionGPList
	do while(associated(fnGP))

!	  if(fnGP%IsProdFunctionGP) then
!	    write(*,*) 'TotalRevCoeffName is ', fnGP%TotalRevCoeff%name
!	  else if(fnGP%IsDemandFunctionGP) then
!	    write(*,*) 'TotalIncomeCoeff is ', fnGP%TotalIncomeCoeff%name
!	  end if

	  fnGP => fnGP%next
	end do

!	write(*,*)' End WriteFunctionGPList'

	end subroutine WriteFunctionGPList

!!	Allocate a new functionGP pointer with default values

	subroutine newfunctionGP(fnGP, fnMPSGE)
! Input
	type (function), pointer ::fnMPSGE
! Output
	type (functionGP), pointer :: fnGP

	nullify(fnGP)
	allocate(fnGP)
	fnGP%functionMPSGE => fnMPSGE
	return
	end subroutine newfunctionGP

	subroutine ConvertToGP

	type (functionGP), pointer :: fnGP
	type (function), pointer :: fnMPSGE
	type (field), pointer :: fnf, NetputIODEField
	type (symbol), pointer :: ThisSymbol
	type (GPCoeffVarReal) :: ImpliedVariable, ImpliedLevelsVariable
	integer :: kpass, NumNetputGPDone, NumEndowmentGPDone
	type (netputGP), pointer :: CurrentNetputGP  
	type (taxGP), pointer :: CurrentTaxGP  

!  Does declarations for $prod
!  Based on Tom Rutherford's routine "declarations"
	type (record), pointer :: r
	type (nest), pointer :: n
	type (netputlist), pointer :: e
	type (nest), pointer :: inputs, outputs
	type (domain), pointer :: InputsSubDomain
! for $prod
	type (GPCoeffVarReal) :: InputPrice, InputQuantity 
	type (GPCoeffVarReal) :: OutputPrice, OutputQuantity 
	type (GPCoeffVarReal) :: ThisProdVar, ThisTopVar
	type (GPCoeffVarReal) :: ThisElasticity
	type (symbol), pointer :: ElasticitySymbol
! for $demand
	type (endowmentlist), pointer :: CurrentEndowmentList
	type (endowmentGP), pointer :: CurrentEndowmentGP
	type (GPCoeffVarReal) :: DemandPrice, DemandQuantity 
	type (GPCoeffVarReal) :: EndowmentPrice, EndowmentQuantity
	type (GPCoeffVarReal) :: ThisDemandVar
	logical :: DFieldExists, EFieldExists, Found 

	integer :: CurrentTaxNumber
	integer :: fnindex = 0
	character (len=1024) :: rec
	logical :: DoThis

	character (len=512) :: StatementStart, OutputString, OutputStringAlls
	character (len=512) :: TempString1

10015   format(1x,a)     ! KP
10017   format(2x,3a)

	write(*,*)' Beginning routine ConvertToGP'

! Whether to translate to a mixed or a linear TAB file
!  Only one of these will be true
!	logical :: MixedTABFile, LinearTABFile
	MixedTABFile = .true.
	LinearTABFile = .false.

! Common declarations
!	write(*,*)'Calling StartTABGP'
	call StartTABGP
!	write(*,*)'Ended StartTABGP'

! Initialize Preline message strings. [Leave the PostLine ones saying are converting to GP.]
	call InitializeMessages('pre',0)

! Run through the sectors, commodities and consumers declaring these as variables
	do kpass=1,3
! Do sectors when KPASS=1
! Do commidities when KPASS=2
! Do consumers when KPASS=3
	  if(kpass.eq.1) then
	    write(iotabgp, '(/,A)') '! $sectors: !'
	    write(iotabgp, '(A)') '!  [The variables below are %-changes in the activity level of these sectors.] !'
	    Messages%PreLine(1) = 'Processing $sectors'
	  else if(kpass.eq.2) then
	    write(iotabgp, '(/,A)') '! $commodities: !'
	    write(iotabgp, '(A)') '!  [These are market prices] !'
	    Messages%PreLine(1) = 'Processing $commodities'
	  else if(kpass.eq.3) then
	    write(iotabgp, '(/,A)') '! $consumers: !'
	    write(iotabgp, '(A)') '!  [The variables below are %-changes in the income of these consumers.] !'
! Don't need to add ORIG_LEVEL for $consumers since later will declare levels
!  variable (total value of income for the relevant $consumer) which is associated with each $consumer   
	    Messages%PreLine(1) = 'Processing $consumers'
	  end if
!  Re-initialize 'pre' messages past number 1
	  call InitializeMessages('pre',1)
 
	  ThisSymbol => idlist    ! idlist contains all symbols
	  do while(associated(ThisSymbol))
	    DoThis = .false.
	    if(kpass.eq.1) DoThis = ThisSymbol%type.eq.'sector'
	    if(kpass.eq.2) DoThis = ThisSymbol%type.eq.'commodity'
	    if(kpass.eq.3) DoThis = ThisSymbol%type.eq.'consumer'
! Sectors, Commodities and Consumers
	    if(DoThis) then

	      Messages%PreLine(2) = 'Processing '//trim(ThisSymbol%name)
	      call InitializeMessages('pre',2)

! Declare this as a variable - use ImpliedVariable to hold info
	      call InitializeGPCoeffVarReal('v', ImpliedVariable)
	      ImpliedVariable%name = ThisSymbol%name
	      ImpliedVariable%arg => ThisSymbol%arg
	      ImpliedVariable%labelinfo = trim(ThisSymbol%labelinfo)//' (%-change)'

	      if(kpass.eq.1) then

! Orig_Level for $sectors. Find this in a $prod: function
!  [The orig_level value for these is the sum of the values in the q: fields, summed over all o: lines.
!    Indeed, the orig_level is coeffficient R_ (TotalRevCoeff).] 

		fnGP => FunctionGPList
		Found = .false.
		do while(associated(fnGP).and. (.not.Found))
		  if(fnGP%IsProdFunctionGP) then
		    if(EqualStringsIgnoreCase(ThisSymbol%name, fnGP%ActivityLevelVar%name)) then
		      Found = .true.
		      ImpliedVariable%OrigLevel = fnGP%TotalRevCoeff%name
		    end if
		  end if
		  fnGP => fnGP%next
		end do ! do while(associated(fnGP).and. (.not.Found))
	      end if ! if(kpass.eq.1) then

! Omit Orig_Level now that also declare levels variables - instead add Formula(Initial)s for
!  the levels variables - see below
! NO	      if(kpass.eq.2) then
!		  Assume all prices are 1 initially
! NO	        ImpliedVariable%OrigLevel = '1.0'
! NO	      end if

	      call DeclareGPCoeffVar(ImpliedVariable)
! Also add this name to GPNameList
	      call AddToGPNameList( ImpliedVariable%name, 'u', 'v', ImpliedVariable%arg)

	      if(kpass.eq.2) then

! Check that name is not too long
		if(LEN_TRIM(ThisSymbol%name)+2.ge.llcs) then
		  Call ErrorStop('Sorry, names of $consumers are limited to '// &
				trim(PosintToString(llcs-2))//' characters.')		
		end if ! if(LEN_TRIM(ThisSymbol%name)+2.ge.llcs) then

! Also declare levels values of prices for commodities
	        call InitializeGPCoeffVarReal('c', ImpliedLevelsVariable)
	        ImpliedLevelsVariable%name = trim(ThisSymbol%name)//'_L'
	        ImpliedLevelsVariable%arg => ThisSymbol%arg
	        ImpliedLevelsVariable%labelinfo = trim(ThisSymbol%labelinfo)//' (levels value)'
		call DeclareGPLevelsVar(ImpliedLevelsVariable, ImpliedVariable)
! Also add this name to GPNameList
	        call AddToGPNameList( ImpliedLevelsVariable%name, 'u', 'c', ImpliedLevelsVariable%arg)

! And give Formula(Initial) for its original value
		StatementStart = 'Formula (Initial)'
		call ArgToAllsGP(ImpliedLevelsVariable%arg, OutputStringAlls)
! Write this line
		OutputString = trim(StatementStart)//' '//trim(OutputStringAlls)
	        call WriteOneLineToTAB(OutputString)
! LHS, "=" and 1.0
		call CoeffVarRealToGP(ImpliedLevelsVariable, TempString1)
		OutputString = trim(TempString1)//' = '//'1.0 ;' 
! Write this line
	        call WriteOneLineToTAB('   '//OutputString)

	      end if ! if(kpass.eq.2) then

	    end if ! if(DoThis) then

	    ThisSymbol => ThisSymbol%next
	  end do
	end do ! kpass

!	write(*,*)'Ended sectors etc'	

	fnGP => FunctionGPList
	do while (associated(fnGP))
	  fnMPSGE => fnGP%functionMPSGE
	  fnf => fnMPSGE%firstrecord%firstfield
	  if ((fnf%label.eq.'$prod').or.(fnf%label.eq.'$demand')) then

	    Messages%PreLine(1) = 'Processing function '//trim(fnf%label)//':'//trim(fnf%data)
!  Re-initialize 'pre' messages past number 1
	    call InitializeMessages('pre',1)
!	    write(*,*) trim(Messages%PreLine(1))

! Put out comment above the declarations for this function (just echos the MPS text)
!	    write(*,*)'Calling OutputMPSFunctionAsGPComment'
	    call OutputMPSFunctionAsGPComment(fnMPSGE)
!	    write(*,*)'Called OutputMPSFunctionAsGPComment'

	  end if

	  if (fnf%label.eq.'$prod') then

! Declare any elasticities associated with this function
!	    write(*,*)'Starting S elasticity'
	    if(fnGP%SElastPresent) then
	      if(.not.fnGP%SElastIsConstant) then
	        call DeclareGPCoeffVar(fnGP%SElasticity)
	      end if
	    end if
!	    write(*,*)'Starting T elasticity'
	    if(fnGP%TElastPresent) then
	      if(.not.fnGP%TElastIsConstant) then
	        call DeclareGPCoeffVar(fnGP%TElasticity)
	      end if
	    end if

!	    write(*,*)'Done elasticities'

! Declarations for:
!	  type (GPCoeffVarReal) :: TotalRevCoeff !! Marginal revenue as coefficient (relates to outputs)
!	  type (GPCoeffVarReal) :: MargRevVar !! Marginal revenue as variable (relates to outputs)
!	  type (GPCoeffVarReal) :: TotalCostCoeff !! Marginal cost as coefficient (relates to inputs)
!	  type (GPCoeffVarReal) :: MargCostVar !! Marginal cost as variable (relates to inputs)
	    call DeclareGPCoeffVar(fnGP%TotalRevCoeff)
	    call DeclareGPCoeffVar(fnGP%TotalCostCoeff)
	    call DeclareGPCoeffVar(fnGP%MargRevVar)
	    call DeclareGPCoeffVar(fnGP%MargCostVar)

!	    write(*,*)'Declared to point 1'

! Declarations for each i: field and each o: field
!		 Set CurrentNetputGP to point to first i: or o: field associated with
!		   this function
	    CurrentNetputGP => fnGP%FirstNetputGP
	    NumNetputGPDone = 0
	    do while(NumNetputGPDone.lt.fnGP%NumIFields+fnGP%NumOFields)

	      if(.not.(associated(CurrentNetputGP))) then
		call InternalErrorStop( &
		  'CurrentNetputGP is not associated. NumNetputGPDone while in routine ConvertToGP')
	      end if


	      NetputIODEField => CurrentNetputGP%NetputMPSGE%f   ! This is the :, o:, d: or e: field at start of this netput
	      Messages%PreLine(2) = 'Processing line beginning '// &
	      				trim(NetputIODEField%label)//':'//trim(NetputIODEField%data)
!  Re-initialize 'pre' messages past number 2
	      call InitializeMessages('pre',2)
!	      write(*,*) trim(Messages%PreLine(2))

!	      if(CurrentNetputGP%QFieldIsExpression) then
!		write(*,*) trim(Messages%PreLine(2)), 'QFieldIsExpression is true'
!	      else
!		write(*,*) trim(Messages%PreLine(2)), 'QFieldIsExpression is false'
!	      end if 


! Do declarations for:
!	  type (GPCoeffVarReal) :: IOQFieldQuantityVar  !! %-change in input or output quantity
!	  type (GPCoeffVarReal) :: IOPriceVar    !! %-change in price
!	  type (GPCoeffVarReal) :: IOQFieldValueCoeff  !! value (level) of input or output
!	  type (GPCoeffVarReal) :: IOQFieldQuantityCoeffOrRealMPSGE  !! MPSGE quantity (initial level) of input or output

! Do not declare IOPriceVar since it should have been declared under $commodities
! NO	      call DeclareGPCoeffVar(CurrentNetputGP%IOPriceVar)

! If writing mixed TAB file, IOQFieldValueCoeff is a levels variable if are any
!   tax terms on this line
! Now always levels if MixedTABFile
!WAS	      if(MixedTABFile.and.CurrentNetputGP%NumTaxTypes.gt.0) then
	      if(MixedTABFile) then
	        call DeclareGPLevelsVar(CurrentNetputGP%IOQFieldValueCoeff, CurrentNetputGP%IOQFieldValueVar)
	      else
	        call DeclareGPCoeffVar(CurrentNetputGP%IOQFieldValueCoeff)
	      end if

! Next seems to be wrong
! NO. IOQFieldQuantityVar is not a new declaration if there is no subdomain
! NO.	      if(associated(CurrentNetputGP%NetputMPSGE%sd)) then

! Is the ORIG_LEVEL in the next variable correct?
	        write(iotabgp, '(a)') ' ! ?? Is the ORIG_LEVEL in next variable correct ?? !'

! IOQFieldQuantityVar is now always associated with IOQuantityCoeff which is a levels variable
! WAS	        call DeclareGPCoeffVar(CurrentNetputGP%IOQFieldQuantityVar)
	      if(MixedTABFile) then
	        call DeclareGPLevelsVar(CurrentNetputGP%IOQFieldQuantityCoeff, CurrentNetputGP%IOQFieldQuantityVar)
	      else
	        call DeclareGPCoeffVar(CurrentNetputGP%IOQFieldQuantityCoeff)
	      end if

! NO.	      else
! Should be able to add ORIG_LEVEL to the quantity associated with the $prod.
! NO.	        write(iotabgp, '(a)') ' ! ?? Should associate ORIG_LEVEL with the activity variable in this $prod ?? !'
! NO.	      end if

! IOQFieldQuantityCoeffOrRealMPSGE - not set up if CurrentNetputGP%QFieldIsExpression is true
	      if(.not.CurrentNetputGP%QFieldIsExpression) then
!  Nothing to declare if this is a real
	       if(CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%type.eq.'c') then
! Only declare it if it is different from IOQFieldValueCoeff (the one updated in GEMPACK)
	        if(trim(CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%name).ne. &
	               trim(CurrentNetputGP%IOQFieldQuantityCoeff%name)) then
	          call DeclareGPCoeffVar(CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE)
	        end if
	       end if
	      end if

	      if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Run through the taxes
		CurrentTaxGP => CurrentNetputGP%FirstTaxGP
		CurrentTaxNumber = 0
		do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)

	          Messages%PreLine(3) = 'Processing tax with agent '//'a:'// &
	      				trim(CurrentTaxGP%taxlistMPSGE%a%data)
!  Re-initialize 'pre' messages past number 3
	          call InitializeMessages('pre',3)
!	          write(*,*) trim(Messages%PreLine(2))

! Do declarations for:
!	  type (GPCoeffVarReal) :: TaxPowerCoeff  !! power of tax
!	  type (GPCoeffVarReal) :: TaxRateCoeff  !! rate of tax (number like 0.13)
!	  type (GPCoeffVarReal) :: TaxRevenueCoeff  !! tax revenue
!	  type (GPCoeffVarReal) :: TaxPowerVar  !! %-change in power of tax
!	  type (GPCoeffVarReal) :: TaxRateChangeVar  !! change in tax rate (number like 0.05)
!	  type (GPCoeffVarReal) :: TaxRevenueChangeVar  !! change in tax revenue
! Next will be a real (eg t:0.2) or a Coefficient (eg t:t_p(o))
!	  type (GPCoeffVarReal) :: TaxRateCoeffOrRealMPSGE  !! MPSGE value (initial level) of tax rate

		  if(MixedTABFile) then
! Levels variables
	            call DeclareGPLevelsVar(CurrentTaxGP%TaxPowerCoeff, CurrentTaxGP%TaxPowerVar)
	            call DeclareGPLevelsVar(CurrentTaxGP%TaxRateCoeff, CurrentTaxGP%TaxRateChangeVar)
	            call DeclareGPLevelsVar(CurrentTaxGP%TaxRevenueCoeff, CurrentTaxGP%TaxRevenueChangeVar)
		  else
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxPowerCoeff)
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxRateCoeff)
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxRevenueCoeff)
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxPowerVar)
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxRateChangeVar)
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxRevenueChangeVar)
		  end if

! TaxRateCoeffOrRealMPSGE
!  Nothing to declare if this is a real
	          if(CurrentTaxGP%TaxRateCoeffOrRealMPSGE%type.eq.'c') then
	            call DeclareGPCoeffVar(CurrentTaxGP%TaxRateCoeffOrRealMPSGE)
		  end if

		  CurrentTaxNumber = CurrentTaxNumber + 1
		  CurrentTaxGP => CurrentTaxGP%next
		end do ! while

	        if(CurrentNetputGP%NumTaxTypes.ge.2) then

! Next only needed if are 2 or more tax types for this netput
!	  type (GPCoeffVarReal) :: TotalTaxPowerCoeff  !! power of tax (all tax types)
!	  type (GPCoeffVarReal) :: TotalTaxPowerVar  !! %-change in power of tax (all tax types)

		  if(MixedTABFile) then
! Levels variable
	            call DeclareGPLevelsVar(CurrentNetputGP%TotalTaxPowerCoeff, CurrentNetputGP%TotalTaxPowerVar)
		  else
	            call DeclareGPCoeffVar(CurrentNetputGP%TotalTaxPowerCoeff)
	            call DeclareGPCoeffVar(CurrentNetputGP%TotalTaxPowerVar)
		  end if

!	         End "if(CurrentNetputGP%NumTaxTypes.ge.2) then"
	        end if

!	  integer :: NumTaxTypes
! Next only needed if are one or more taxing agents
!	  type (GPCoeffVarReal) :: IOPriceIncTaxesVar    !! %-change in price including taxes
! Next only needed for outputs when are one or more taxing agents. 
!	  type (GPCoeffVarReal) :: IOQFieldValueVar    !! %-change in value of input or output

	        write(iotabgp, '(a)') ' ! ?? Is the ORIG_LEVEL in next variable correct ?? !'
	        call DeclareGPCoeffVar(CurrentNetputGP%IOPriceIncTaxesVar)

! If writing a mixed TAB file, this is part of a levels variable already declared
			
		if(LinearTABFile.and.CurrentNetputGP%NetputMPSGE%f%label.eq.'o') then
	          write(iotabgp, '(a)') ' ! ?? Is the ORIG_LEVEL in next variable correct ?? !'
	          call DeclareGPCoeffVar(CurrentNetputGP%IOQFieldValueVar)
		end if

!	       End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
	      end if
	      
! Increment NumNetputGPDone
	      NumNetputGPDone = NumNetputGPDone + 1

! Go on to the next netput for this function
	      CurrentNetputGP => CurrentNetputGP%next

!	     End "do while(NumNetputGPDone.lt.fnGP%NumIFields+fnGP%NumOFields)"
	    end do

!	   End "if (fnf%label.eq.'$prod') then"
	  else if(fnf%label.eq.'$demand') then

! Declarations for:
!	  type (GPCoeffVarReal) :: IncomeCoeff !! total income for this agent

	    if(MixedTABFile) then
! Levels variable
	      call DeclareGPLevelsVar(fnGP%TotalIncomeCoeff, FnGP%IncomeVar)
	    else
	      call DeclareGPCoeffVar(FnGP%TotalIncomeCoeff)
	    end if
! Don't declare IncomeVar since it has been declared in $consumers block


	    if(fnGP%NumEFields.gt.0) then

! Declarations for each e: field
!		 Set CurrentEndowmentGP to point to first e: field associated with
!		   this function
	      CurrentEndowmentGP => fnGP%FirstEndowmentGP
	      NumEndowmentGPDone = 0
	      do while(NumEndowmentGPDone.lt.fnGP%NumEFields)

	        Messages%PreLine(2) = 'Processing line beginning '// &
	      				'e:'//trim(CurrentEndowmentGP%endowmentMPSGE%f%data)
!  Re-initialize 'pre' messages past number 2
	        call InitializeMessages('pre',2)
!	        write(*,*) trim(Messages%PreLine(2))

!	        if(CurrentEndowmentGP%QFieldIsExpression) then
!		  write(*,*) trim(Messages%PreLine(2)), 'QFieldIsExpression is true'
!	        else
!		  write(*,*) trim(Messages%PreLine(2)), 'QFieldIsExpression is false'
!	        end if 
! Declarations for:
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeff  !! quantity (level) of each endowment q: field  
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityVar  !! change in quantity of each endowment (q: field) 
!	  type (GPCoeffVarReal) :: EndowQFieldValueCoeff  !! value (level) of endowment
!	  type (GPCoeffVarReal) :: EndowQFieldValueVar  !! %-change in value (level) of endowment

! Next will be a real (eg q:200) or a Coefficient (eg q:totfac)
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeffOrRealMPSGE  !! MPSGE quantity (initial level) of endowment

	        if(MixedTABFile) then
! Levels variable
	          call DeclareGPLevelsVar(CurrentEndowmentGP%EndowQFieldQuantityCoeff, &
						CurrentEndowmentGP%EndowQFieldQuantityVar)
		else
		  call DeclareGPCoeffVar(CurrentEndowmentGP%EndowQFieldQuantityVar)
		end if

	        if(MixedTABFile) then
! Levels variable
	          call DeclareGPLevelsVar(CurrentEndowmentGP%EndowQFieldValueCoeff, CurrentEndowmentGP%EndowQFieldValueVar)
		else
	          call DeclareGPCoeffVar(CurrentEndowmentGP%EndowQFieldValueCoeff)
		end if

!  EndowQFieldQuantityCoeffOrRealMPSGE is not set up if CurrentEndowmentGP%QFieldIsExpression is true
		if(.not.CurrentEndowmentGP%QFieldIsExpression) then
!  Nothing to declare if this is a real
	         if(CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%type.eq.'c') then

! Only declare it if it is different from EndowQFieldQunatityCoeff (the one updated in GEMPACK)
! If this e: line has an r: field in it, need to compare EndowQFieldQuantityCoeffOrRealMPSGE with EndowUnrationedQuantityCoeff
		  DoThis = .false.
		  if(CurrentEndowmentGP%HasRField) then
	           DoThis = .not.EqualStringsIgnoreCase(CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%name, &
	               CurrentEndowmentGP%EndowUnrationedQuantityCoeff%name)
		    write(*,*)'CurrentEndowmentGP%EndowUnrationedQuantityCoeff is ', &
			CurrentEndowmentGP%EndowUnrationedQuantityCoeff%name
		  else
	           DoThis = .not.EqualStringsIgnoreCase(CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%name, &
	               CurrentEndowmentGP%EndowQFieldQuantityCoeff%name)
		    write(*,*)'CurrentEndowmentGP%EndowQFieldQuantityCoeff is ', &
			CurrentEndowmentGP%EndowQFieldQuantityCoeff%name
		  end if
		  if(DoThis) then
		    write(*,*)'DoThis is true. CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE is ', &
			CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%name
		  else
		    write(*,*)'DoThis is false. CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE is ', &
			CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%name
		  end if 


	          if(DoThis) then
	            call DeclareGPCoeffVar(CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE)
	          end if

	         end if
	        end if

! EndowRFieldCoeff  !! level of r: variable
		if(CurrentEndowmentGP%HasRField) then
	          if(MixedTABFile) then
! Levels variable
	            call DeclareGPLevelsVar(CurrentEndowmentGP%EndowRFieldCoeff, CurrentEndowmentGP%EndowRFieldVar)
		  else
	            call DeclareGPCoeffVar(CurrentEndowmentGP%EndowRFieldCoeff)
		  end if

	          if(MixedTABFile) then
! Levels variable
	            call DeclareGPLevelsVar(CurrentEndowmentGP%EndowUnRationedQuantityCoeff, &
							CurrentEndowmentGP%EndowUnRationedQuantityVar)
		  else
	            call DeclareGPCoeffVar(CurrentEndowmentGP%EndowUnRationedQuantityCoeff)
		  end if

		end if ! if(CurrentEndowmentGP%HasRField) then

! Increment NumEndowmentGPDone
	        NumEndowmentGPDone = NumEndowmentGPDone + 1

	        CurrentEndowmentGP => CurrentEndowmentGP%next
	      end do ! while

!	     End "if(fnGP%NumEFields.gt.0) then"
	    end if


!	   End "else if(fnf%label.eq.'$demand') then"
	  end if


	  if(fnf%label.eq.'$prod'.or.fnf%label.eq.'$demand') then
! Do other statements associated with this function
	    call DoRestTABOneFunctionGP(fnGP)
	  end if


	  fnGP => fnGP%next
	end do

! Write the market clearing equations
	Messages%PreLine(1) = 'Doing the market clearing equations'
!  Re-initialize 'pre' messages past number 1
	call InitializeMessages('pre',1)
!	write(*,*) trim(Messages%PreLine(1))
	call DoMarketClearingEquations

! Process the constraints
	Messages%PreLine(1) = 'Doing the constraints'
!  Re-initialize 'pre' messages past number 1
	call InitializeMessages('pre',1)
	write(*,*) trim(Messages%PreLine(1))
	call DoConstraints

! Blank line at end of TAB file
	call WriteOneLineToTAB(' ')


! Write top of CMF file
	Messages%PreLine(1) = 'Doing the Command file'
!  Re-initialize 'pre' messages past number 1
	call InitializeMessages('pre',1)
	write(*,*) trim(Messages%PreLine(1))
	call WriteCMFFile

! Close and keep the output files
	close(unit=iotabgp, status='keep')
	close(unit=iocmfgp, status='keep')

	return
	end subroutine ConvertToGP

	subroutine DoRestTABOneFunctionGP(ThisFunctionGP)
! Input
	type (functionGP), pointer :: ThisFunctionGP
! Local
	type (function), pointer :: ThisFunctionMPSGE
	integer :: kpass, NumNetputGPDone, kcoeff, CurrentTaxNumber, NumEndowmentGPDone
	type (netputGP), pointer :: CurrentNetputGP
	type (taxGP), pointer :: CurrentTaxGP
	type (endowmentGP), pointer :: CurrentEndowmentGP
	character(len=512) OutputString, OutputStringAlls, OutputString1, OutputString2
	character(len=512) OutputStringArgs, CoeffVarName, DeclarationStart
	integer :: LenStart, LenSoFar, LenAlls
	logical :: RHSDone

10005	format(a)

	ThisFunctionMPSGE => ThisFunctionGP%functionMPSGE

	do kpass=1,4
! Do Reads when KPASS=1
! Do Formulas when KPASS=2
! Do Equations when KPASS=3
! Do Updates when KPASS=4

	  if (ThisFunctionMPSGE%firstrecord%firstfield%label.eq.'$prod') then

	    if(kpass.eq.1) then
	      call DoReadsOneProdFunction(ThisFunctionGP)
	    end if

	    CurrentNetputGP => ThisFunctionGP%FirstNetputGP
	    NumNetputGPDone = 0
	    do while(NumNetputGPDone.lt.ThisFunctionGP%NumIFields+ThisFunctionGP%NumOFields)

! [All reads for this function have already been done]
	      if(kpass.eq.2) then
! In the mixedTABFile case, DoFormulasOneIOLine may write "formula & equation" statements
	        call DoFormulasOneIOLine(CurrentNetputGP)
	      else if(kpass.eq.3) then
	        call DoEquationsOneIOLine(ThisFunctionGP, CurrentNetputGP)
	      else if(kpass.eq.4) then
	        call DoUpdatesOneIOLine(CurrentNetputGP)
	      end if


! Increment NumNetputGPDone
	      NumNetputGPDone = NumNetputGPDone + 1

! Go on to the next netputGP for this function
	      CurrentNetputGP => CurrentNetputGP%next

	    end do ! while

	    if(kpass.eq.2) then
! [06/04/04] Need to delay these until KPASS=3 so that values on RHS are set up
! Formulas for Total revenue and cost
!LATER	      call DoCostRevenueFormulasOneProdFunction(ThisFunctionGP)
	    else if(kpass.eq.3) then
! Formulas for Total revenue and cost
	      call DoCostRevenueFormulasOneProdFunction(ThisFunctionGP)
! Equations for total revenue and cost, and zero profit equation
	      call DoEquationsOneProdFunction(ThisFunctionGP)
	    end if

!	   End "if (ThisFunctionMPSGE%firstrecord%firstfield%label.eq.'$prod') then"
	  else if (ThisFunctionMPSGE%firstrecord%firstfield%label.eq.'$demand') then

	    if(ThisFunctionGP%NumEFields.gt.0) then
	      CurrentEndowmentGP => ThisFunctionGP%FirstEndowmentGP
	      NumEndowmentGPDone = 0
	      do while(NumEndowmentGPDone.lt.ThisFunctionGP%NumEFields)

	        if(kpass.eq.1) then
!  Read value in q: field
	          call DoReadsOneEndowmentLine(CurrentEndowmentGP)
		else if(kpass.eq.2) then
	          call DoFormulasOneEndowmentLine(CurrentEndowmentGP)
		else if(kpass.eq.3) then
!  Linear equation for %-change in value of q: field
	          call DoEquationsOneEndowmentLine(CurrentEndowmentGP)
	        end if

! Increment NumEndowmentGPDone
	        NumEndowmentGPDone = NumEndowmentGPDone + 1

! Go on to the next endowmentGP for this function
	        CurrentEndowmentGP => CurrentEndowmentGP%next
	      end do ! while

!	     End "if(ThisFunctionGP%NumEIFields.gt.0) then"
	    end if

	    if(kpass.eq.3) then
!		Levels equation for total income - must be done after DoEquationsOneEndowmentLine
	      call DoEquationsOneDemandFunction(ThisFunctionGP)
	    end if

!	   End "else if (ThisFunctionMPSGE%firstrecord%firstfield%label.eq.'$demand') then"
	  end if

	end do ! kpass=1,4

	return
	end subroutine DoRestTABOneFunctionGP

	recursive subroutine WriteOneLineToTAB(ThisLine) 
!! Routine WriteOneLineToTAB: Write one line, breaking it up into lines of at most 80 character if necessary
! Input
	character (len=*), intent(in) :: ThisLine
! Local
	integer :: n1
	character (len=512) :: TempLine

10005	format(a)

	n1 = LEN_TRIM(ThisLine)
	if(n1.le.80) then
	  write(iotabgp, 10005) trim(ThisLine)
	else
! Try to break the line at sensible places.
! First try "="
	  n1 = INDEX(ThisLine,'=')
	  if(n1.eq.0) then
! Next try "+"
	    n1=INDEX(ThisLine,'+')
	  end if
	  if(n1.eq.0) then
! Next try "-"
	    n1=INDEX(ThisLine,'-')
	  end if
	  if(n1.eq.0) then
! Next try "*"
	    n1=INDEX(ThisLine,'*')
	  end if
	  if(n1.gt.0) then
! Write as 2 lines
! First line ends at "=","+" or "-".
	    call WriteOneLineToTAB(ThisLine(1:n1))
	    call WriteOneLineToTAB(ThisLine(n1+1:))
	  else
! Try to break before things like SUM
	    TempLine = ThisLine
	    call lower(TempLine)
	    n1=INDEX(TempLine,'sum')
	    if(n1.gt.1) then
! Write as 2 lines [Note that the first part may be all blanks]
	      TempLine = ThisLine(1:n1-1)
	      if(TempLine.ne.' ') then
	        call WriteOneLineToTAB(TempLine(1:n1-1))
	      end if
	      call WriteOneLineToTAB(ThisLine(n1:))
	    else
	      call InternalErrorStop( &
	        'WriteOneLineToTAB: Line is too long. Length is:'//PosintToString(n1) // &
		  '\n Line is:'// &
	            '\n'//trim(ThisLine) )
	    end if
	  end if
	end if

	return
	end subroutine WriteOneLineToTAB 

	subroutine DoUpdatesOneIOLine(CurrentNetputGP)
! Input
	type (netputGP), pointer :: CurrentNetputGP
! Local
	integer :: kcoeff, CurrentTaxNumber
	type (taxGP), pointer :: CurrentTaxGP
	character(len=512) OutputString, OutputStringAlls, OutputString1, OutputString2
	character(len=512) OutputStringArgs, CoeffVarName, DeclarationStart
	integer :: LenStart, LenSoFar, LenAlls
	logical :: RHSDone

10005	format(a)

! Taxes
	if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Run through the taxes
		CurrentTaxGP => CurrentNetputGP%FirstTaxGP
		CurrentTaxNumber = 0
		do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)

! Update TaxRateCoeff  !! rate of tax (number like 0.13)
! Update TaxRevenueCoeff  !! tax revenue

! in MixedTABFile case, never need to write these updates
		  if(LinearTABFile) then
  
		    do kcoeff=1,2
! Do TaxRateCoeff when KCOEFF=1
! Do TaxRevenueCoeff when KCOEFF=2

		      DeclarationStart = 'Update (Change)'
		      LenStart = LEN_TRIM(DeclarationStart)

		      OutputString=trim(DeclarationStart)//' '
		      LenSoFar = LenStart+1

		      call ArgToAllsGP(CurrentTaxGP%TaxRateCoeff%arg, OutputStringAlls)
		      OutputStringAlls = trim(OutputStringAlls)
		      LenAlls = LEN_TRIM(OutputStringAlls)
		      if(OutputStringAlls.ne.'') then
! Output Start now if ncessary
	  	        if(LenStart.gt.25) then
			  call WriteOneLineToTAB(OutputString)
	    		  OutputString = ' '
	  	        end if
	  	        OutputString = trim(OutputString)//' '//OutputStringAlls
		      end if
		      LenSoFar = LEN_TRIM(OutputString)

! Output a line now if ncessary
		      if(LenSoFar.gt.40) then
			call WriteOneLineToTAB(OutputString)
	  	        OutputString = ' '
		      end if

! The LHS name
		      if(kcoeff.eq.1) then
		        CoeffVarName = CurrentTaxGP%TaxRateCoeff%name
		      else if(kcoeff.eq.2) then
		        CoeffVarName = CurrentTaxGP%TaxRevenueCoeff%name
		      end if
		      OutputString = trim(OutputString)//' '//trim(CoeffVarName)
! The arguments
		      OutputStringArgs = ''
!
		      if(OutputStringAlls.ne.'') then

		        call ArgToGP(CurrentTaxGP%TaxRateCoeff%arg, OutputStringArgs)
! Add opening ( and closing )
		        OutputStringArgs = '('//trim(OutputStringArgs)//')'
! Add to OutputString
		        OutputString = trim(OutputString)//OutputStringArgs
		      end if
! '='
		      OutputString = trim(OutputString)//' = '

! Output a line now if ncessary
		      if(LenSoFar.gt.40) then
			call WriteOneLineToTAB(OutputString)
	  	        OutputString = ' '
		      end if

! RHS name
		      if(kcoeff.eq.1) then
		        CoeffVarName = CurrentTaxGP%TaxRateChangeVar%name
		      else if(kcoeff.eq.2) then
		        CoeffVarName = CurrentTaxGP%TaxRevenueChangeVar%name
		      end if
		      OutputString = trim(OutputString)//' '//trim(CoeffVarName)

! RHS arguments (same as on LHS in this case) and ";"
		      OutputString = trim(OutputString)//trim(OutputStringArgs)//' ;'		

! Complete the statement
		      call WriteOneLineToTAB(OutputString)

		    end do ! kcoeff=1,2

!		   End "if(LinearTABFile) then"
		  end if

		  CurrentTaxNumber = CurrentTaxNumber + 1
		  CurrentTaxGP => CurrentTaxGP%next
		end do ! while

!	       End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
	end if

! Update IOQFieldValueCoeff
		
! This Update is not needed in MixedTABFile case if are any taxes on this line
! [06/04/04] Nor is it ever needed in MixedTABFile case
!WAS	if(LinearTABFile .or. CurrentNetputGP%NumTaxTypes.eq.0) then
	if(LinearTABFile) then
 
		DeclarationStart = 'Update'
		LenStart = LEN_TRIM(DeclarationStart)

		OutputString=trim(DeclarationStart)//' '
		LenSoFar = LenStart+1

		call ArgToAllsGP(CurrentNetputGP%IOQFieldValueCoeff%arg, OutputStringAlls)
		OutputStringAlls = trim(OutputStringAlls)
		LenAlls = LEN_TRIM(OutputStringAlls)
		if(OutputStringAlls.ne.'') then
! Output Start now if ncessary
	  	        if(LenStart.gt.25) then
			  call WriteOneLineToTAB(OutputString)
	    		  OutputString = ' '
	  	        end if
	  	        OutputString = trim(OutputString)//' '//OutputStringAlls
		end if
		LenSoFar = LEN_TRIM(OutputString)

! Output a line now if ncessary
		if(LenSoFar.gt.40) then
			call WriteOneLineToTAB(OutputString)
	  	        OutputString = ' '
		end if

! The LHS name
		CoeffVarName = CurrentNetputGP%IOQFieldValueCoeff%name
		OutputString = trim(OutputString)//' '//trim(CoeffVarName)
! The arguments
		OutputStringArgs = ''
!
		if(OutputStringAlls.ne.'') then

		  call ArgToGP(CurrentNetputGP%IOQFieldValueCoeff%arg, OutputStringArgs)
! Add opening ( and closing )
		  OutputStringArgs = '('//trim(OutputStringArgs)//')'
! Add to OutputString
		  OutputString = trim(OutputString)//OutputStringArgs
		end if
! '='
		OutputString = trim(OutputString)//' = '

! Output a line now if ncessary
		if(LenSoFar.gt.40) then
			call WriteOneLineToTAB(OutputString)
	  	        OutputString = ' '
		end if

! RHS
		RHSDone = .false.
		if(CurrentNetputGP%NumTaxTypes.ge.1) then
		  if(CurrentNetputGP%netputMPSGE%f%label.eq.'o') then
! Use IOQFieldValueVar to do the update
		    CoeffVarName = CurrentNetputGP%IOQFieldValueVar%name
		    OutputString = trim(OutputString)//' '//trim(CoeffVarName)
! RHS arguments (same as on LHS in this case) and ";"
		    OutputString = trim(OutputString)//trim(OutputStringArgs)//' ;'		
! Complete the statement
		    call WriteOneLineToTAB(OutputString)
		    RHSDone = .true .
		  end if
		end if

		if(.not.RHSDone) then

! Use product update, price variable is IOPriceVar and quantity variable is IOQFieldQuantityVar 
		  call CoeffVarRealToGP(CurrentNetputGP%IOPriceVar, OutputString1) 
		  call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityVar, OutputString2)
		  OutputString = trim(OutputString)//' '//trim(OutputString1)// &
				' * '//trim(OutputString2)//' ;' 
! Complete the statement
		  call WriteOneLineToTAB(OutputString)
		  RHSDone = .true .

!		 End "if(.not.RHSDone) then"
		end if

	end if ! if(LinearTABFile) then

	end subroutine DoUpdatesOneIOLine

	subroutine DoOneReadIfNotAlreadyDone(CoeffNameToRead)
! input
	character (len=*), intent(in) :: CoeffNameToRead
! Local
	character(len=512) OutputString
	character (len=llhd) HeadString
	logical :: found

10005	format(a)
 
! Has this already been read?
	call IsNameInGPReadList(CoeffNameToRead, Found)
	if(.not.found) then
! Find header to read this from
	  call GetReadHeader(CoeffNameToRead, HeadString)
	  OutputString = 'Read '//trim(CoeffNameToRead)// &
		  ' from File Input Header "'//HeadString//'" ;'
	  call WriteOneLineToTAB(OutputString)
! Add this read to GPReadList
	  call AddToGPReadList(CoeffNameToRead)
	end if

	end subroutine DoOneReadIfNotAlreadyDone

	subroutine DoReadsOneProdFunction(ThisFunctionGP)
! Input
	type (functionGP), pointer :: ThisFunctionGP
! Local
	integer :: NumNetputGPDone
	type (netputGP), pointer :: CurrentNetputGP
	character(len=512) OutputString, ElastName
	character (len=4) HeadString
	logical :: found, DoThis
	type (GPCoeffVarReal) :: ThisElast
	integer :: kpass

10005	format(a)
 
	CurrentNetputGP => ThisFunctionGP%FirstNetputGP
	NumNetputGPDone = 0
	do while(NumNetputGPDone.lt.ThisFunctionGP%NumIFields+ThisFunctionGP%NumOFields)

	  call DoReadsOneIOLine(CurrentNetputGP)

! Increment NumNetputGPDone
	  NumNetputGPDone = NumNetputGPDone + 1

! Go on to the next netput for this function
	  CurrentNetputGP => CurrentNetputGP%next

	end do ! while

! Read the elasticities if they are not reals
	do kpass=1,2

! SElast when KPASS=1
! TElast when KPASS=2
	  if(kpass.eq.1) then
	    DoThis = ThisFunctionGP%SElastPresent
	  else if(kpass.eq.2) then
	    DoThis = ThisFunctionGP%TElastPresent
	  end if
	  if(DoThis) then
	    if(kpass.eq.1) then
	      ThisElast = ThisFunctionGP%SElasticity
	    else if(kpass.eq.2) then
	      ThisElast = ThisFunctionGP%TElasticity
	    end if
! Nothing to read if this is a real
	    if(ThisElast%type.eq.'c') then
	      call DoOneReadIfNotAlreadyDone(ThisElast%name)
	    end if	   
	  end if
	end do ! kpass

	return
	end subroutine DoReadsOneProdFunction

	subroutine DoReadsOneIOLine(CurrentNetputGP)
! Input
	type (netputGP), pointer :: CurrentNetputGP
! Local
	integer :: CurrentTaxNumber
	type (taxGP), pointer :: CurrentTaxGP

10005	format(a)

! Taxes
	if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Run through the taxes
	  CurrentTaxGP => CurrentNetputGP%FirstTaxGP
	  CurrentTaxNumber = 0
	  do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)

! Read TaxRateCoeffOrRealMPSGE if it is not a real number
	    if(CurrentTaxGP%TaxRateCoeffOrRealMPSGE%type.eq.'c') then
	      call DoOneReadIfNotAlreadyDone(CurrentTaxGP%TaxRateCoeffOrRealMPSGE%name)
	    end if

	    CurrentTaxNumber = CurrentTaxNumber + 1
	    CurrentTaxGP => CurrentTaxGP%next

	  end do ! while

!	       End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
	end if

! Read IOQFieldQuantityCoeffOrRealMPSGE if it is not a real number and if it is set up
	if(.not.CurrentNetputGP%QFieldIsExpression) then
	 if(CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%type.eq.'c') then
	  call DoOneReadIfNotAlreadyDone(CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%name)
	 end if
	end if
		
	return
	end subroutine DoReadsOneIOLine

	subroutine DoReadsOneEndowmentLine(CurrentEndowmentGP)
! Input
	type (endowmentGP), pointer :: CurrentEndowmentGP

10005	format(a)

! Read EndowQFieldQuantityCoeffOrRealMPSGE if it is not a real number and it is set up
	if(.not.CurrentEndowmentGP%QFieldIsExpression) then
	 if(CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%type.eq.'c') then
	  call DoOneReadIfNotAlreadyDone(CurrentEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%name)
	 end if
	end if

! Read EndowRFieldCoeff if HasRField is true
	if(CurrentEndowmentGP%HasRField) then
	  call DoOneReadIfNotAlreadyDone(CurrentEndowmentGP%EndowRFieldCoeff%name)
	end if
		
	return
	end subroutine DoReadsOneEndowmentLine

	subroutine DoFormulasOneEndowmentLine(ThisEndowmentGP)
! Input
	type(endowmentGP), pointer :: ThisEndowmentGP
! Local
	logical :: DoThis
	type (GPCoeffVarReal) :: LHSCoeff  
	character(len=512) OutputString, TempString, TempString2

! Formula(Initial) for EndowQFieldQuantityCoeff or EndowUnrationedQuantityCoeff if q: field is an
!  expression or a number

	DoThis = .false.
	if(ThisEndowmentGP%QFieldIsExpression) then
	  DoThis = .true.
	else
	  DoThis = .not.(ThisEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE%type.eq.'c')
	end if

	if(DoThis) then

!	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeff  !! quantity (level) of each endowment q: field  
!	  type (GPCoeffVarReal) :: EndowUnrationedQuantityCoeff  !! unrationed quantity of endowment (level)

! LHS coefficient is EndowQFieldQuantityCoeff if no r: field, otherwise is EndowUnrationedQuantityCoeff
	  if(ThisEndowmentGP%HasRField) then
	    LHSCoeff = ThisEndowmentGP%EndowUnrationedQuantityCoeff
	  else
	    LHSCoeff = ThisEndowmentGP%EndowQFieldQuantityCoeff
	  end if

! Formula keyword and quantifiers
	  call ArgToAllsGP(LHSCoeff%arg, TempString)
	  OutputString = 'Formula (Initial) '//trim(TempString)
	  call WriteOneLineToTAB(OutputString)
	
! The LHS and =
	  call CoeffVarRealToGP(LHSCoeff, TempString)
	  OutputString = trim(TempString)//' = '
	  call WriteOneLineToTAB(OutputString)

! RHS and ";"
! RHS is QFieldExpressionString when QFieldIsExpression is true
	  if(ThisEndowmentGP%QFieldIsExpression) then
	    OutputString = trim(ThisEndowmentGP%QFieldExpressionString)//' ;'
	  else
! The RHS: IOQFieldQuantityCoeffOrRealMPSGE
	    call CoeffVarRealToGP(ThisEndowmentGP%EndowQFieldQuantityCoeffOrRealMPSGE, TempString)
	    OutputString = trim(OutputString)//' '//trim(TempString)//' ;'
	  end if
! Write this out
	  call WriteOneLineToTAB('  '//OutputString)

	end if ! if(DoThis) then

	return
	end subroutine DoFormulasOneEndowmentLine

	subroutine DoFormulasOneIOLine(CurrentNetputGP)
! In the mixedTABFile case, DoFormulasOneIOLine may write "formula & equation" statements
! Input
	type (netputGP), pointer :: CurrentNetputGP
! Local
	integer :: kcoeff, CurrentTaxNumber
	type (taxGP), pointer :: CurrentTaxGP
	character(len=512) OutputString, TempString, TempString2
	character (len=4) HeadString
	logical :: found, IsOutputLine, DoThis
	character :: PlusMinusChar  ! either '+' or '-', depending on IsOutputLine
	character (len=80) :: StatementStart

10005	format(a)

! Need Formula (Initial) for value of q: field if either QFieldIsExpression, or if
!   IOQFieldQuantityCoeffOrRealMPSGE is different from IOQFieldQuantityCoeff).
	DoThis = CurrentNetputGP%QFieldIsExpression .or. &
	    .not.EqualStringsIgnoreCase(CurrentNetputGP%IOQFieldQuantityCoeff%name, &
			 CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE%name)

	if(DoThis) then

! Formula(Initial) for IOQFieldQuantityCoeff - equal to IOQFieldValueCoeffMPSGE (if they are not the same) 
!	  type (GPCoeffVarReal) :: IOQFieldQuantityCoeff  !! quantity (level) of input or output
! Next will be a real (eg q:200) or a Coefficient (eg q:totfac)
!	  type (GPCoeffVarReal) :: IOQFieldQuantityCoeffOrRealMPSGE  !! MPSGE quantity (initial level) of input or output

! Formula keyword and quantifiers
	  call ArgToAllsGP(CurrentNetputGP%IOQFieldQuantityCoeff%arg, TempString)
	  OutputString = 'Formula (Initial) '//trim(TempString)
	  call WriteOneLineToTAB(OutputString)
! The LHS: IOQFieldQuantityCoeff
	  call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityCoeff, TempString)
	  OutputString = trim(TempString)//' = '
	  if(LEN_TRIM(OutputString).ge.40) then
! Write this out
	    call WriteOneLineToTAB('  '//OutputString)
	    OutputString = ''
	  end if
! RHS is QFieldExpressionString when QFieldIsExpression is true
	  if(CurrentNetputGP%QFieldIsExpression) then
	    OutputString = trim(OutputString)//' '//trim(CurrentNetputGP%QFieldExpressionString)//' ;'
	  else
! The RHS: IOQFieldQuantityCoeffOrRealMPSGE
	    call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityCoeffOrRealMPSGE, TempString)
	    OutputString = trim(OutputString)//' '//trim(TempString)//' ;'
	  end if
! Write this out
	  call WriteOneLineToTAB('  '//OutputString)

	end if ! if(DoThis) then


! Formula & Equation saying that IOQFieldValueVar = IOPriceVar * IOQFieldQuantityVar
	StatementStart = 'Formula & Equation'
	OutputString = trim(StatementStart)//' E_'//CurrentNetputGP%IOQFieldValueVar%name
! Write this out
	call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
	call MarkAsEndogenous(CurrentNetputGP%IOQFieldValueVar)
! Quantifiers and LHS and '='
	call ArgToAllsGP(CurrentNetputGP%IOQFieldValueCoeff%arg, TempString)
	call CoeffVarRealToGP(CurrentNetputGP%IOQFieldValueCoeff, TempString2)
	OutputString = trim(TempString)//' '//trim(TempString2)//' ='
! Write this out
	call WriteOneLineToTAB('  '//OutputString)

! RHS
	call CoeffVarRealToGP(CurrentNetputGP%IOPriceCoeff, TempString)
	call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityCoeff, TempString2)
	OutputString = trim(TempString)//' * '//trim(TempString2)//' ;'

! Write this out
	call WriteOneLineToTAB('    '//OutputString)

! IsOutputLine and PlusMinusChar
	IsOutputLine = CurrentNetputGP%netputMPSGE%f%label.eq.'o'
	if(IsOutputLine) then
	  PlusMinusChar = '-'
	else
	  PlusMinusChar = '+'
	end if

! Taxes
	if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Run through the taxes
	  CurrentTaxGP => CurrentNetputGP%FirstTaxGP
	  CurrentTaxNumber = 0
	  do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)

! Do formulas associated with this tax type
!	  type (GPCoeffVarReal) :: TaxPowerCoeff  !! power of tax
!	  type (GPCoeffVarReal) :: TaxRateCoeff  !! rate of tax (number like 0.13)
!	  type (GPCoeffVarReal) :: TaxRevenueCoeff  !! tax revenue

! Formula(Initial) for TaxRateCoeff
! Formula keyword and quantifiers
	    call ArgToAllsGP(CurrentTaxGP%TaxRateCoeff%arg, TempString)
	    OutputString = 'Formula (Initial) '//trim(TempString)
	    call WriteOneLineToTAB(OutputString)
! The LHS
	    call CoeffVarRealToGP(CurrentTaxGP%TaxRateCoeff, TempString)
	    OutputString = trim(TempString)//' = '
	    if(LEN_TRIM(OutputString).gt.40) then
! Write this out
	      call WriteOneLineToTAB('  '//OutputString)
	      OutputString = ''
	    end if
! The RHS which is TaxRateCoeffOrRealMPSGE 
	    call CoeffVarRealToGP(CurrentTaxGP%TaxRateCoeffOrRealMPSGE, TempString)
	    OutputString = trim(OutputString)//' '//trim(TempString)//' ;'
	    call WriteOneLineToTAB('  '//OutputString)

! Formula or "Formula & Equation" for TaxPowerCoeff
! Formula keyword and quantifiers
	    call ArgToAllsGP(CurrentTaxGP%TaxPowerCoeff%arg, TempString)
! Set up StatementStart
	    if(MixedTABFile) then
! Add the equation name here
	      StatementStart = 'Formula & Equation '//'E_'//CurrentTaxGP%TaxPowerVar%name
	    else
	      StatementStart = 'Formula'
	    end if
	    OutputString = trim(StatementStart)//' '//trim(TempString)
	    call WriteOneLineToTAB(OutputString)
! Mark the RATE (not the POWER) as endogenous
	    call MarkAsEndogenous(CurrentTaxGP%TaxRateChangeVar)
! The LHS
	    call CoeffVarRealToGP(CurrentTaxGP%TaxPowerCoeff, TempString)
	    OutputString = trim(TempString)//' = '
	    if(LEN_TRIM(OutputString).gt.40) then
! Write this out
	      call WriteOneLineToTAB('  '//OutputString)
	      OutputString = ''
	    end if
! The RHS which is 1 - TaxRateCoeff [Output line] or 1 + TaxRateCoeff [Input line] 
	    call CoeffVarRealToGP(CurrentTaxGP%TaxRateCoeff, TempString)
	    OutputString = trim(OutputString)//' 1 '//PlusMinusChar//' '//trim(TempString)//' ;'
	    call WriteOneLineToTAB('  '//OutputString)

! Formula (Initial) or "Formula & Equation" for TaxRevenueCoeff
! Formula keyword and quantifiers
	    call ArgToAllsGP(CurrentTaxGP%TaxRevenueCoeff%arg, TempString)
! Set up StatementStart
	    if(MixedTABFile) then
! Add the equation name here
	      StatementStart = 'Formula & Equation '//'E_'//CurrentTaxGP%TaxRevenueChangeVar%name
	    else
	      StatementStart = 'Formula (Initial)'
	    end if
	    OutputString = trim(StatementStart)//' '//trim(TempString)
	    call WriteOneLineToTAB(OutputString)
! Mark this as endogenous
	    call MarkAsEndogenous(CurrentTaxGP%TaxRevenueChangeVar)
! The LHS
	    call CoeffVarRealToGP(CurrentTaxGP%TaxRevenueCoeff, TempString)
	    OutputString = trim(TempString)//' = '
	    if(LEN_TRIM(OutputString).gt.30) then
! Write this out
	      call WriteOneLineToTAB('  '//OutputString)
	      OutputString = ''
	    end if
! The RHS which is TaxRateCoeff * IOQFieldValueCoeff 
	    call CoeffVarRealToGP(CurrentTaxGP%TaxRateCoeff, TempString)
	    call CoeffVarRealToGP(CurrentNetputGP%IOQFieldValueCoeff, TempString2)
	    OutputString = trim(OutputString)//' '//trim(TempString)//' * '//trim(TempString2)//' ;'
	    call WriteOneLineToTAB('  '//OutputString)

	    CurrentTaxNumber = CurrentTaxNumber + 1
	    CurrentTaxGP => CurrentTaxGP%next

!	   End "do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)"
	  end do ! while

	  if(CurrentNetputGP%NumTaxTypes.ge.2) then

! Formula [(Always)] or "Formula & Equation" for TotalTaxPowerCoeff
!	  type (GPCoeffVarReal) :: TotalTaxPowerCoeff  !! power of tax (all tax types)

! Formula keyword and quantifiers
	    call ArgToAllsGP(CurrentNetputGP%TotalTaxPowerCoeff%arg, TempString)
! Set up StatementStart
	    if(MixedTABFile.and.CurrentNetputGP%NumTaxTypes.ge.2) then
! Add the equation name here
	      StatementStart = 'Formula & Equation '//'E_'//CurrentNetputGP%TotalTaxPowerVar%name
	    else
	      StatementStart = 'Formula'
	    end if
	    OutputString = trim(StatementStart)//' '//trim(TempString)
	    call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
	    call MarkAsEndogenous(CurrentNetputGP%TotalTaxPowerVar)
! The LHS
	    call CoeffVarRealToGP(CurrentNetputGP%TotalTaxPowerCoeff, TempString)
	    OutputString = trim(TempString)//' = '
! Write this out
	    call WriteOneLineToTAB('  '//OutputString)

! The RHS which is 1 several [+/- TaxRateCoeff] terms
	    OutputString = '1'  
! Run through the taxes
	    CurrentTaxGP => CurrentNetputGP%FirstTaxGP
	    CurrentTaxNumber = 0
	    do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)

!	  type (GPCoeffVarReal) :: TaxRateCoeff  !! tax rate
	      call CoeffVarRealToGP(CurrentTaxGP%TaxRateCoeff, TempString)
	      OutputString = trim(OutputString)//' '//PlusMinusChar//' '//TempString
	      if(LEN_TRIM(OutputString).gt.40) then
! Write this out
	        call WriteOneLineToTAB('  '//OutputString)
	        OutputString = ''
	      end if

	      CurrentTaxNumber = CurrentTaxNumber + 1
	      CurrentTaxGP => CurrentTaxGP%next

!	     End "do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)"
	    end do ! while

! Semicolon to end the formula
	    OutputString = trim(OutputString)//' ;'
! Write this line to complete the formula
	    call WriteOneLineToTAB('  '//OutputString)

!	   End "if(CurrentNetputGP%NumTaxTypes.ge.2) then"
	  end if

!	       End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
	end if

	return
	end subroutine DoFormulasOneIOLine

	subroutine DoCostRevenueFormulasOneProdFunction(ThisFunctionGP)
! Input
	type (functionGP), pointer :: ThisFunctionGP
! Local
	type (netputGP), pointer :: ThisNetputGP
	type (taxGP), pointer :: ThisTaxGP
	integer :: NumFieldsToDo, NumFieldsDone, NSum, kpass
	type (GPCoeffVarReal) :: ThisGPCoeff, ThisGPCoeff2, LHSGPCoeff
	logical :: UseThis
	type (domain), pointer :: ThisSubDomain
	character(len=512) OutputString, TempString, TempString2, SumString, EndSumString
	integer :: CurrentTaxNumber
	character :: PlusOrMinusChar  ! '-' for outputs, '+' for inputs

	do kpass=1,2

! Outputs when KPASS=1
! Inputs when KPASS=2

!	  type (GPCoeffVarReal) :: TotalRevCoeff !! Marginal revenue as coefficient (relates to outputs)
!	  type (GPCoeffVarReal) :: TotalCostCoeff !! Marginal cost as coefficient (relates to inputs)

	  if(kpass.eq.1) then
	    LHSGPCoeff = ThisFunctionGP%TotalRevCoeff
	    NumFieldsToDo = ThisFunctionGP%NumOFields
	    PlusOrMinusChar = '-'
	  else if(kpass.eq.2) then
	    LHSGPCoeff = ThisFunctionGP%TotalCostCoeff
	    NumFieldsToDo = ThisFunctionGP%NumIFields
	    PlusOrMinusChar = '+'
	  end if

! Formula keyword and quantifiers
	  call ArgToAllsGP(LHSGPCoeff%arg, TempString)
	  OutputString = 'Formula '//trim(TempString)
	  call WriteOneLineToTAB(OutputString)
! The LHS
	  call CoeffVarRealToGP(LHSGPCoeff, TempString)
	  OutputString = trim(TempString)//' = '
! Write this out
	  call WriteOneLineToTAB('  '//OutputString)

! RHS.
! Must add tax revenue (for inputs) or subtract tax revenue (for outputs)

! Go through all the relevant (input or output) fields
	  NumFieldsDone = 0
	  ThisNetputGP => ThisFunctionGP%FirstNetputGP
	  do while(NumFieldsDone.lt.NumFieldsToDo)

! Only use this if is of correct (input or output) type
	    if(kpass.eq.1) then
	      UseThis = ThisNetputGP%netputMPSGE%f%label.eq.'o'
	    else if(kpass.eq.2) then
	      UseThis = ThisNetputGP%netputMPSGE%f%label.eq.'i'
	    end if

	    if(UseThis) then

! Revenue or Cost for this line
! If is a subdomain, need to sum over it

	      ThisSubdomain => ThisNetputGP%netputMPSGE%sd
	      call WriteSumsGP(ThisSubDomain, SumString, NSum)
	      call CoeffVarRealToGP(ThisNetputGP%IOQFieldValueCoeff, TempString)
! End the sums
	      call EndSumsGP(NSum, EndSumString)

	      if(ThisNetputGP%NumTaxTypes.eq.0) then

! No tax terms to consider
! Assemble and write it out for this output line
	        OutputString = trim(SumString)//' '//trim(TempString)//' '//trim(EndSumString)
		if(NumFieldsDone.le.NumFieldsToDo-2) then
!		     "+" at end before next term
		    OutputString = trim(OutputString)//' +'
		else
! semicolon to end the formula
		    OutputString = trim(OutputString)//' ;'
		end if

! Write this out
	        call WriteOneLineToTAB('    '//OutputString)

	      else if(ThisNetputGP%NumTaxTypes.gt.0) then

! Must consider the tax terms

! First output SUM and IOQFieldValueCoeff
	        OutputString = trim(SumString)//' '//trim(TempString)
	        call WriteOneLineToTAB('    '//OutputString)
		
! Run through the taxes
	        ThisTaxGP => ThisNetputGP%FirstTaxGP
	        CurrentTaxNumber = 0
	        do while(CurrentTaxNumber.lt.ThisNetputGP%NumTaxTypes)

	          call CoeffVarRealToGP(ThisTaxGP%TaxRevenueCoeff, TempString)
	          OutputString = PlusOrMinusChar//' '//TempString
		  if(CurrentTaxNumber.eq.ThisNetputGP%NumTaxTypes-1) then

!			This is the last term. Add the ENDSUMs here.
		    OutputString = trim (OutputString)//' '//EndSumString

		    if(NumFieldsDone.le.NumFieldsToDo-2) then
!		       "+" at end before next term
		      OutputString = trim(OutputString)//' +'
		    else
!		       ";" to end formula
		      OutputString = trim(OutputString)//' ;'
		    end if

		  end if


! Write this out
	          call WriteOneLineToTAB('        '//OutputString)

	          CurrentTaxNumber = CurrentTaxNumber + 1
	          ThisTaxGP => ThisTaxGP%next

!	         End "do while(CurrentTaxNumber.lt.ThisNetputGP%NumTaxTypes)"
	        end do ! while

!	       End "if(ThisNetputGP%NumTaxTypes.eq.0) then..else.."
	      end if

	      NumFieldsDone = NumFieldsDone + 1

!	     End "if(UseThis) then"
	    end if

	    ThisNetputGP => ThisNetputGP%next

!	   End "do while(NumFieldsDone.le.NumFieldsToDo)"
	  end do

	end do ! kpass=1,2

	return
	end subroutine DoCostRevenueFormulasOneProdFunction

	subroutine DoEquationsOneIOLine(CurrentFunctionGP, CurrentNetputGP)
! Input
	type (functionGP), pointer :: CurrentFunctionGP
	type (netputGP), pointer :: CurrentNetputGP
! Local
	character(len=512) OutputString, StatementStart, TempString, TempString2
	logical :: IsOLine, IsILine
	character PlusOrMinusChar

10005	format(a)

! Linear CES/CET equation associated with this i: or o: line
!  LHS is IOQFieldQuantityVar
	if(MixedTABFile) then
	  StatementStart = 'Equation (Linear)'
	else
	  StatementStart = 'Equation'
	end if
	OutputString = trim(StatementStart)//' E_'//CurrentNetputGP%IOQFieldQuantityVar%name
! Write this out
	call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
	call MarkAsEndogenous(CurrentNetputGP%IOQFieldQuantityVar)

! Quantifiers and LHS and '='
	call ArgToAllsGP(CurrentNetputGP%IOQFieldQuantityVar%arg, TempString)
	call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityVar, TempString2)
	OutputString = trim(TempString)//' '//trim(TempString2)//' ='
! Write this out
	call WriteOneLineToTAB('  '//OutputString)

! RHS which is:
!  ActivityLevelVar + TElasticity * [<var>  - MargRevVar ]  for an o: line 
!  ActivityLevelVar - SElasticity * [<var>  - MargCostVar ]  for an i: line
! Above <var> is IOPriceVar if no taxes on this line, or is IOPriceIncTaxesVar if are some taxes  
	IsOLine = CurrentNetputGP%netputMPSGE%f%label.eq.'o'
	IsILine = CurrentNetputGP%netputMPSGE%f%label.eq.'i'

	if(IsOLine) then
	  PlusOrMinusChar = '+'
	else
	  PlusOrMinusChar = '-'
	end if

	call CoeffVarRealToGP(CurrentFunctionGP%ActivityLevelVar, OutputString)
	OutputString = trim(OutputString)//' '//PlusOrMinusChar
	if(IsOLine) then
	  if(CurrentFunctionGP%TElastPresent) then
	    call CoeffVarRealToGP(CurrentFunctionGP%TElasticity, TempString)
	  else
	    TempString = '0'
	  end if
	else
	  if(CurrentFunctionGP%SElastPresent) then
	    call CoeffVarRealToGP(CurrentFunctionGP%SElasticity, TempString)
	  else
	    TempString = '0'
	  end if
	end if
	OutputString = trim(OutputString)//' '//trim(TempString)//' * ['
	if(CurrentNetputGP%NumTaxTypes.gt.0) then
	  call CoeffVarRealToGP(CurrentNetputGP%IOPriceIncTaxesVar, TempString)
	else
	  call CoeffVarRealToGP(CurrentNetputGP%IOPriceVar, TempString)
	end if
	OutputString = trim(OutputString)//trim(TempString)//' -'
  
	if(IsOLine) then
	  call CoeffVarRealToGP(CurrentFunctionGP%MargRevVar, TempString)
	else
	  call CoeffVarRealToGP(CurrentFunctionGP%MargCostVar, TempString)
	end if
	OutputString = trim(OutputString)//' '//trim(TempString)//'] ;'
! Write this out
	call WriteOneLineToTAB('    '//OutputString)

! Taxes
! Need to do this for OLine or ILine if are taxes
! [06/04/04] This is now always done, and in the levels (see above)
!NOWLEV	if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Equation for %-change in value of q: field. [This is only for o: lines]
! LHS is IOQFieldValueVar
! RHS is always Market Price (IOPriceVar) + IOQFieldQuantityVar. This is true for o: and i: fields

!NOWLEV	  if(MixedTABFile) then
!NOWLEV	    StatementStart = 'Equation (Linear)'
!NOWLEV	  else
!NOWLEV	    StatementStart = 'Equation'
!NOWLEV	  end if
!NOWLEV	  OutputString = trim(StatementStart)//' E_'//CurrentNetputGP%IOQFieldValueVar%name
! Write this out
!NOWLEV	  call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
!NOWLEV	  call MarkAsEndogenous(CurrentNetputGP%IOQFieldValueVar)

! Quantifiers and LHS and '='
!NOWLEV	  call ArgToAllsGP(CurrentNetputGP%IOQFieldValueVar%arg, TempString)
!NOWLEV	  call CoeffVarRealToGP(CurrentNetputGP%IOQFieldValueVar, TempString2)
!NOWLEV	  OutputString = trim(TempString)//' '//trim(TempString2)//' ='
! Write this out
!NOWLEV	  call WriteOneLineToTAB('  '//OutputString)

! RHS
!NOWLEV	  call CoeffVarRealToGP(CurrentNetputGP%IOPriceVar, TempString)
!NOWLEV	  call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityVar, TempString2)
!NOWLEV	  OutputString = trim(TempString)//' + '//trim(TempString2)//' ;'

! Write this out
!NOWLEV	  call WriteOneLineToTAB('    '//OutputString)

!	 End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
!NOWLEV	end if 
	
	if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Equation connecting price in i: or o: field with price including taxes and power of tax
! LHS is IOPriceIncTaxesVar
! RHS is IOPriceVar + <var> where
!    <var> is TotalTaxPowerVar if are 2 or more tax types, or
!    <var> is TaxPowerVar for the 1 tax if is just one tax type

	  if(MixedTABFile) then
	    StatementStart = 'Equation (Linear)'
	  else
	    StatementStart = 'Equation'
	  end if
	  OutputString = trim(StatementStart)//' E_'//CurrentNetputGP%IOPriceIncTaxesVar%name
! Write this out
	  call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
	call MarkAsEndogenous(CurrentNetputGP%IOPriceIncTaxesVar)

! Quantifiers and LHS and '='
	  call ArgToAllsGP(CurrentNetputGP%IOPriceIncTaxesVar%arg, TempString)
	  call CoeffVarRealToGP(CurrentNetputGP%IOPriceIncTaxesVar, TempString2)
	  OutputString = trim(TempString)//' '//trim(TempString2)//' ='
! Write this out
	  call WriteOneLineToTAB('  '//OutputString)

! RHS
	  call CoeffVarRealToGP(CurrentNetputGP%IOPriceVar, TempString)
	  if(CurrentNetputGP%NumTaxTypes.ge.2) then
	    call CoeffVarRealToGP(CurrentNetputGP%TotalTaxPowerVar, TempString2)
	  else     !if(CurrentNetputGP%NumTaxTypes.eq.1) then
	    call CoeffVarRealToGP(CurrentNetputGP%FirstTaxGP%TaxPowerVar, TempString2)
	  end if
	  OutputString = trim(TempString)//' + '//trim(TempString2)//' ;'
! Write this out
	  call WriteOneLineToTAB('    '//OutputString)

!	 End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
	end if

!------------------------------------------------------------
! *** Next not necessary *** ??
!??	if(CurrentNetputGP%NumTaxTypes.gt.0) then

! Run through the taxes
!??	  CurrentTaxGP => CurrentNetputGP%FirstTaxGP
!??	  CurrentTaxNumber = 0
!??	  do while(CurrentTaxNumber.lt.CurrentNetputGP%NumTaxTypes)


!??	    CurrentTaxNumber = CurrentTaxNumber + 1
!??	    CurrentTaxGP => CurrentTaxGP%next

!??	  end do ! while

!	       End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
!??	end if
! ** End of not necessary? part ***
!------------------------------------------------------------

	end subroutine DoEquationsOneIOLine

	subroutine DoEquationsOneProdFunction(ThisFunctionGP)
! Input
	type (functionGP), pointer :: ThisFunctionGP
! Local
	character(len=512) OutputString, StatementStart, TempString, TempString2
	character(len=512) TempString3, TempString4, TempString5
	character(len=512) PriceVarString, SumString, EndSumString
	character :: PlusOrMinusChar  ! '-' for outputs, '+' for inputs
	type (netputGP), pointer :: ThisNetputGP
	type (taxGP), pointer :: ThisTaxGP
	integer :: CurrentTaxNumber, NumFieldsDone, NumFieldsToDo, kpass, NSum
	logical :: UseThis
	type (domain), pointer :: ThisSubDomain
	character :: ISChar

	do kpass=1,2

! Equation for marginal revenue (output lines) when KPASS=1
! Equation for marginal cost (input lines) when KPASS=2

! This is based on the code which writes the formulas for TotalRevenue and 
!   TotalCost in routine DoCostRevenueFormulasOneProdFunction

!	  type (GPCoeffVarReal) :: TotalRevCoeff !! Marginal revenue as coefficient (relates to outputs)
!	  type (GPCoeffVarReal) :: MargRevVar !! Marginal revenue as variable (relates to outputs)
!	  type (GPCoeffVarReal) :: TotalCostCoeff !! Total cost as coefficient (relates to inputs)
!	  type (GPCoeffVarReal) :: MargCostVar !! Marginal cost as variable (relates to inputs)

	  if(MixedTABFile) then
	    StatementStart = 'Equation (Linear)'
	  else
	    StatementStart = 'Equation'
	  end if
	  if(kpass.eq.1) then
	    OutputString = trim(StatementStart)//' E_'//ThisFunctionGP%MargRevVar%name
! Mark this variable as endogenous
	    call MarkAsEndogenous(ThisFunctionGP%MargRevVar)
	  else if(kpass.eq.2) then
	    OutputString = trim(StatementStart)//' E_'//ThisFunctionGP%MargCostVar%name
! Mark this variable as endogenous
	    call MarkAsEndogenous(ThisFunctionGP%MargCostVar)
	  end if
! Write this out
	  call WriteOneLineToTAB(OutputString)

! Quantifiers and LHS and '='
	  if(kpass.eq.1) then
	    call ArgToAllsGP(ThisFunctionGP%MargRevVar%arg, OutputString)
	    call CoeffVarRealToGP(ThisFunctionGP%TotalRevCoeff, TempString)
	    call CoeffVarRealToGP(ThisFunctionGP%MargRevVar, TempString2)
	  else if(kpass.eq.2) then
	    call ArgToAllsGP(ThisFunctionGP%MargCostVar%arg, OutputString)
	    call CoeffVarRealToGP(ThisFunctionGP%TotalCostCoeff, TempString)
	    call CoeffVarRealToGP(ThisFunctionGP%MargCostVar, TempString2)
	  end if
	  OutputString = trim(OutputString)//' '//trim(TempString)//' * '//trim(TempString2)//' ='
! Write this out
	  call WriteOneLineToTAB('  '//OutputString)


	  if(kpass.eq.1) then
	    NumFieldsToDo = ThisFunctionGP%NumOFields
	    PlusOrMinusChar = '-'
	  else if(kpass.eq.2) then
	    NumFieldsToDo = ThisFunctionGP%NumIFields
	    PlusOrMinusChar = '+'
	  end if

! RHS.
! Must add tax revenue (for inputs) or subtract tax revenue (for outputs)

! Go through all the relevant (input or output) fields
	  NumFieldsDone = 0
	  ThisNetputGP => ThisFunctionGP%FirstNetputGP
	  do while(NumFieldsDone.lt.NumFieldsToDo)

! Only use this if is of correct (input or output) type
	    if(kpass.eq.1) then
	      UseThis = ThisNetputGP%netputMPSGE%f%label.eq.'o'
	    else if(kpass.eq.2) then
	      UseThis = ThisNetputGP%netputMPSGE%f%label.eq.'i'
	    end if

	    if(UseThis) then

! Revenue or Cost for this line
! If is a subdomain, need to sum over it

	      ThisSubdomain => ThisNetputGP%netputMPSGE%sd
	      call WriteSumsGP(ThisSubDomain, SumString, NSum)
	      call CoeffVarRealToGP(ThisNetputGP%IOQFieldValueCoeff, TempString)
! End the sums
	      call EndSumsGP(NSum, EndSumString)

	      if(ThisNetputGP%NumTaxTypes.eq.0) then

! No tax terms to consider
! Relevant price variable is the price in the o: or i: field
	        call CoeffVarRealToGP(ThisNetputGP%IOPriceVar, PriceVarString)
! Assemble and write it out for this output line
	        OutputString = trim(SumString)//' '//trim(TempString)// &
			  ' *'//trim(PriceVarString)//' '//trim(EndSumString)
		if(NumFieldsDone.le.NumFieldsToDo-2) then
!		     "+" at end before next term
		    OutputString = trim(OutputString)//' +'
		else
! semicolon to end the equation
		    OutputString = trim(OutputString)//' ;'
		end if

! Write this out
	        call WriteOneLineToTAB('    '//OutputString)

	      else if(ThisNetputGP%NumTaxTypes.gt.0) then

! Must consider the tax terms

! Relevant price variable is the price including taxes
	        call CoeffVarRealToGP(ThisNetputGP%IOPriceIncTaxesVar, PriceVarString)

! First output SUM, then opening "[" and then the IOQFieldValueCoeff
	        OutputString = trim(SumString)//' ['//trim(TempString)
	        call WriteOneLineToTAB('    '//OutputString)
		
! Run through the taxes
	        ThisTaxGP => ThisNetputGP%FirstTaxGP
	        CurrentTaxNumber = 0
	        do while(CurrentTaxNumber.lt.ThisNetputGP%NumTaxTypes)

	          call CoeffVarRealToGP(ThisTaxGP%TaxRevenueCoeff, TempString)
	          OutputString = PlusOrMinusChar//' '//TempString
		  if(CurrentTaxNumber.eq.ThisNetputGP%NumTaxTypes-1) then

!			This is the last term. Closing "]", then * PriceVar and the ENDSUMs here.
		    OutputString = trim (OutputString)//'] * '//trim(PriceVarString)//' '//EndSumString

		    if(NumFieldsDone.le.NumFieldsToDo-2) then
!		       "+" at end before next term
		      OutputString = trim(OutputString)//' +'
		    else
!		       ";" to end equation
		      OutputString = trim(OutputString)//' ;'
		    end if

		  end if


! Write this out
	          call WriteOneLineToTAB('        '//OutputString)

	          CurrentTaxNumber = CurrentTaxNumber + 1
	          ThisTaxGP => ThisTaxGP%next

!	         End "do while(CurrentTaxNumber.lt.ThisNetputGP%NumTaxTypes)"
	        end do ! while

!	       End "if(ThisNetputGP%NumTaxTypes.eq.0) then..else.."
	      end if

	      NumFieldsDone = NumFieldsDone + 1

!	     End "if(UseThis) then"
	    end if

	    ThisNetputGP => ThisNetputGP%next

!	   End "do while(NumFieldsDone.le.NumFieldsToDo)"
	  end do

	end do ! kpass=1,2

! Zero profit equation
! LHS is TotalRevCoeff * MargRevVar 
! RHS is TotalCostCoeff * MargCostVar

	if(MixedTABFile) then
	    StatementStart = 'Equation (Linear)'
	else
	    StatementStart = 'Equation'
	end if
! This determines the variable after $prod:
	OutputString = trim(StatementStart)//' E_'// &
	   trim(ThisFunctionGP%ActivityLevelVar%name)// &
	   ' # Zero profit condition for $prod:'// &
	     trim(ThisFunctionGP%functionMPSGE%firstrecord%firstfield%identtext)// &
	   ' #'
! Write this out
	call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
	call MarkAsEndogenous(ThisFunctionGP%ActivityLevelVar)

! Quantifiers
	call ArgToAllsGP(ThisFunctionGP%MargRevVar%arg, OutputString)
! LHS
	call CoeffVarRealToGP(ThisFunctionGP%TotalRevCoeff, TempString)
	call CoeffVarRealToGP(ThisFunctionGP%MargRevVar, TempString2)
! RHS
	call CoeffVarRealToGP(ThisFunctionGP%TotalCostCoeff, TempString3)
	call CoeffVarRealToGP(ThisFunctionGP%MargCostVar, TempString4)

	OutputString = trim(OutputString)//' '//trim(TempString)//' * '// &
	  trim(TempString2)//' = '//trim(TempString3)//' * '//trim(TempString4)//' ;'
! Write this out
	call WriteOneLineToTAB('  '//OutputString)

! Do slightly stricter assertion for intial data, less strict at other steps
	do kpass=1,2
! Do Assertion (Inital) when KPASS=1	
! Do Assertion (Always) when KPASS=2
	  if(kpass.eq.1) then
	    ISChar = 'I'
	  else if(kpass.eq.2) then
	    ISChar = 'S'
	  end if	
! Assertion that costs=revenue
	  if(kpass.eq.1) then
	    StatementStart = 'Assertion (Initial)'
	  else if(kpass.eq.2) then
	    StatementStart = 'Assertion (Always)'
	  end if
	  OutputString = trim(StatementStart)//' # Costs = revenue for $prod:'// &
	     trim(ThisFunctionGP%functionMPSGE%firstrecord%firstfield%identtext)// &
	   ' #'
! Write this out
	  call WriteOneLineToTAB(OutputString)
! Quantifiers
	  call ArgToAllsGP(ThisFunctionGP%MargRevVar%arg, OutputString)
! Write this out
	  call WriteOneLineToTAB(OutputString)
! LHS
	  call CoeffVarRealToGP(ThisFunctionGP%TotalRevCoeff, TempString)
! RHS
	  call CoeffVarRealToGP(ThisFunctionGP%TotalCostCoeff, TempString3)

! Cannot expect exact equality. So instead write it as
!  Revenue > 0.9999*Costs and Revenue < 1.0001*Costs
! [Equality]	OutputString = trim(OutputString)//' '//trim(TempString)// &
! [Equality]	  ' = '//trim(TempString3)//' ;'
! Use BALTOLI or BALTOLS
!	  OutputString = trim(OutputString)//' '//trim(TempString)// &
!	    ' GT 0.9999*'//trim(TempString3)//' AND '//trim(TempString)// &
!	    ' LT 1.0001*'//trim(TempString3)//' ;'
	  TempString4 = '['//trim(TempString)// ' GT (1-BALTOL'//ISChar//')*'//trim(TempString3)//'] AND'
!	  write(*,*)'TempString4 = ', trim(TempString4)
	  TempString2 = '   '//trim(TempString4)
	  call WriteOneLineToTAB(TempString2)
	  TempString5 = '['//trim(TempString)// ' LT (1+BALTOL'//ISChar//')*'//trim(TempString3)//'] ;'
!	  write(*,*)'TempString5 = ', trim(TempString5)
	  TempString2 = '   '//trim(TempString5)
	  call WriteOneLineToTAB(TempString2)
	end do ! do kpass=1,2
	
	return
	end subroutine DoEquationsOneProdFunction

	subroutine PrepareLabelInfoForTABFile(ThisLabellingInfo, OutputString)
! Puts "#" around label and truncates so OutputString has at most 79 characters
! Input
	character (len=*), intent(in) :: ThisLabellingInfo
! Output
	character (len=*), intent(out) :: OutputString
! Local
	character (len=73) :: LabelInfoUse
! Write the label
	LabelInfoUse = trim(ThisLabellingInfo)
! LabelInfoUse has at most 73 characters. Fit label info and ";" on one line
!  function adjustl removes any leading blanks and moves them to the end,
!    but it doesn't seem to work
	OutputString = ' # '//trim(adjustl(LabelInfoUse))//' # '

	return
	end subroutine PrepareLabelInfoForTABFile
	
	subroutine DoEquationsOneEndowmentLine(ThisEndowmentGP)
! Input
	type(endowmentGP), pointer :: ThisEndowmentGP
! Local
	character(len=512) OutputString, StatementStart, TempString, TempString2

	write(*,*)'Beginning DoEquationsOneEndowmentLine'

	if(ThisEndowmentGP%HasRField) then
 
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityCoeff  !! quantity (level) of each endowment q: field  
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityVar  !! change in quantity of each endowment (q: field) 
!	  type (GPCoeffVarReal) :: EndowUnrationedQuantityCoeff  !! unrationed quantity of endowment (level)
!	  type (GPCoeffVarReal) :: EndowUnrationedQuantityVar  !! change in unrationed quantity of endowment

! If is rationing in this line, levels equation saying that 
!   EndowQFieldQuantityCoeff = EndowUnrationedQuantityCoeff * EndowRFieldCoeff (r: field levels variable)

	  if(MixedTABFile) then
	    StatementStart = 'Formula & Equation'
	  else
	    StatementStart = 'Equation'
	  end if
	  OutputString = trim(StatementStart)//' E_'//ThisEndowmentGP%EndowQFieldQuantityVar%name
! Write this out
	  call WriteOneLineToTAB(OutputString)
! Write label
	  call PrepareLabelInfoForTABFile( &
	    'Actual supply connected to unrationed quantity and rationing variable', OutputString)
	  call WriteOneLineToTAB(' '//OutputString)
! Mark this variable as endogenous
	  call MarkAsEndogenous(ThisEndowmentGP%EndowQFieldQuantityVar)

! LHS EndowQFieldQuantityCoeff 
! Quantifiers and LHS and '='
	  call ArgToAllsGP(ThisEndowmentGP%EndowQFieldQuantityCoeff%arg, OutputString)
	  call CoeffVarRealToGP(ThisEndowmentGP%EndowQFieldQuantityCoeff, TempString)
	  OutputString = trim(OutputString)//' '//trim(TempString)//' ='
! Write this out
	  call WriteOneLineToTAB('   '//OutputString)

! RHS. EndowUnrationedQuantityCoeff * EndowRFieldCoeff 
	  call CoeffVarRealToGP(ThisEndowmentGP%EndowUnrationedQuantityCoeff, TempString)
	  call CoeffVarRealToGP(ThisEndowmentGP%EndowRFieldCoeff, TempString2)
	  OutputString = trim(TempString)//' * '//trim(TempString2)//' ;'
! Write this out
	  call WriteOneLineToTAB('          '//OutputString)

	end if ! if(ThisEndowmentGP%HasRField) then

! Formula & Equation saying that EndowQFieldValueCoeff = EndowPriceCoeff * EndowQFieldQuantityCoeff
	StatementStart = 'Formula & Equation'
	OutputString = trim(StatementStart)//' E_'//ThisEndowmentGP%EndowQFieldValueVar%name
! Write this out
	call WriteOneLineToTAB(OutputString)
! Mark this variable as endogenous
	call MarkAsEndogenous(ThisEndowmentGP%EndowQFieldValueVar)
! Quantifiers and LHS and '='
	call ArgToAllsGP(ThisEndowmentGP%EndowQFieldValueCoeff%arg, TempString)
	call CoeffVarRealToGP(ThisEndowmentGP%EndowQFieldValueCoeff, TempString2)
	OutputString = trim(TempString)//' '//trim(TempString2)//' ='
! Write this out
	call WriteOneLineToTAB('  '//OutputString)

! RHS
	call CoeffVarRealToGP(ThisEndowmentGP%EndowPriceCoeff, TempString)
	call CoeffVarRealToGP(ThisEndowmentGP%EndowQFieldQuantityCoeff, TempString2)
	OutputString = trim(TempString)//' * '//trim(TempString2)//' ;'

! Write this out
	call WriteOneLineToTAB('    '//OutputString)



! {06/04/04] This is now always a levels equation - see earlier in this routine
! This writes the equation linking EndowQFieldValueVar with the associated
!   price EndowPriceVar and quantity EndowQFieldQuantityVar
!	  type (GPCoeffVarReal) :: EndowQFieldQuantityVar  !! %-change in quantity of each endowment  
!	  type (GPCoeffVarReal) :: EndowPriceVar    !! %-change in price (The price after e: at start of line)
!	  type (GPCoeffVarReal) :: EndowQFieldValueVar  !! %-change in value (level) of endowment

!LEV	if(MixedTABFile) then
!LEV	  StatementStart = 'Equation (Linear)'
!LEV	else
!LEV	  StatementStart = 'Equation'
!LEV	end if
!LEV	OutputString = trim(StatementStart)//' E_'//ThisEndowmentGP%EndowQFieldValueVar%name
! Write this out
!LEV	call WriteOneLineToTAB(OutputString)
! Write label
!LEV	call PrepareLabelInfoForTABFile( &
!LEV	  'Link '//trim(ThisEndowmentGP%EndowQFieldValueCoeff%name)//' [%-change is '// &
!LEV		trim(ThisEndowmentGP%EndowQFieldValueVar%name)//'] to its price and quantity', &
!LEV		OutputString)
!LEV	call WriteOneLineToTAB(' '//OutputString)
! Mark this variable as endogenous
!LEV	call MarkAsEndogenous(ThisEndowmentGP%EndowQFieldValueVar)

! Quantifiers and LHS and '='
!LEV	call ArgToAllsGP(ThisEndowmentGP%EndowQFieldValueVar%arg, OutputString)
!LEV	call CoeffVarRealToGP(ThisEndowmentGP%EndowQFieldValueVar, TempString)
!LEV	OutputString = trim(OutputString)//' '//trim(TempString)//' ='
! Write this out
!LEV	call WriteOneLineToTAB('   '//OutputString)

! RHS. EndowPriceVar + EndowQFieldQuantityVar
!LEV	call CoeffVarRealToGP(ThisEndowmentGP%EndowPriceVar, TempString)
!LEV	call CoeffVarRealToGP(ThisEndowmentGP%EndowQFieldQuantityVar, TempString2)
!LEV	OutputString = trim(TempString)//' + '//trim(TempString2)//' ;'
! Write this out
!LEV	call WriteOneLineToTAB('          '//OutputString)

	write(*,*)'Ended DoEquationsOneEndowmentLine'

	return
	end subroutine DoEquationsOneEndowmentLine

	subroutine DoEquationsOneDemandFunction(ThisDemandFunctionGP)
! Input
	type(functionGP), pointer :: ThisDemandFunctionGP
! Local
	type (endowmentGP), pointer :: CurrentEndowmentGP
	type (domain), pointer :: ThisSubDomain
	type (functionGP), pointer :: CurrentFunctionGP
	type (netputGP), pointer :: CurrentNetputGP
	type (taxGP), pointer :: CurrentTaxGP
	type (taxlist), pointer :: CurrentTaxMPSGE
	integer :: NumEndowmentGPDone, NSum, NumTaxDone, NumNetputGPDone
	character(len=512) OutputString, StatementStart, TempString
	character(len=512) SumString, EndSumString, ThisConsumer
	integer :: NumTermsSoFar

	ThisConsumer = ThisDemandFunctionGP%functionMPSGE%firstrecord%firstfield%identtext

	NumTermsSoFar = 0

! Just one equation (a levels one if are doing mixed tab file).
! This adds up total income - value of endowments plus value of taxes
! This only done here properly for the MixedTABFile case.

	if(MixedTABFile) then
	  StatementStart = 'Formula & Equation'
	  OutputString = trim(StatementStart)//' E_'//ThisDemandFunctionGP%IncomeVar%name
! Mark this variable as endogenous
	  call MarkAsEndogenous(ThisDemandFunctionGP%IncomeVar)
	else
	    OutputString = 'Formula'
	end if
! Write this out
	call WriteOneLineToTAB(OutputString)
! Write label
	call PrepareLabelInfoForTABFile( &
		ThisDemandFunctionGP%TotalIncomeCoeff%labelinfo, OutputString)
	call WriteOneLineToTAB(' '//OutputString)

! Quantifiers and LHS and '='
	call ArgToAllsGP(ThisDemandFunctionGP%TotalIncomeCoeff%arg, OutputString)
	call CoeffVarRealToGP(ThisDemandFunctionGP%TotalIncomeCoeff, TempString)
	OutputString = trim(OutputString)//' '//trim(TempString)//' ='
! Write this out
	call WriteOneLineToTAB('  '//OutputString)

! RHS
! First go through the endowments (if any)
	if(ThisDemandFunctionGP%NumEFields.gt.0) then

	  CurrentEndowmentGP => ThisDemandFunctionGP%FirstEndowmentGP
	  NumEndowmentGPDone = 0
	  do while(NumEndowmentGPDone.lt.ThisDemandFunctionGP%NumEFields)

! Add value of these endowments
	    if(NumTermsSoFar.gt.0) then
	      StatementStart = '+'
	    else
	      StatementStart = ''
	    end if

! may involve SUM
	    ThisSubdomain => CurrentEndowmentGP%endowmentMPSGE%sd
	    call WriteSumsGP(ThisSubDomain, SumString, NSum)
	    call CoeffVarRealToGP(CurrentEndowmentGP%EndowQFieldValueCoeff, TempString)
! End the sums
	    call EndSumsGP(NSum, EndSumString)
	    OutputString = trim(StatementStart)//' '//trim(SumString)// &
		' '//trim(TempString)//' '//EndSumString 
! Write this out
	    call WriteOneLineToTAB('    '//OutputString)
! Increment NumTermsSoFar
	    NumTermsSoFar = NumTermsSoFar + 1

	    CurrentEndowmentGP => CurrentEndowmentGP%next
	    NumEndowmentGPDone = NumEndowmentGPDone + 1
	  end do ! while

!	 End "if(ThisDemandFunction%NumEFields.gt.0) then"
	end if

! Now add any tax revenue accruing to this agent
! Run through the $prod functionGPs
	CurrentFunctionGP => FunctionGPList
	do while(associated(CurrentFunctionGP))
	  if(CurrentFunctionGP%IsProdFunctionGP) then

! Run through the o: and i: lines for this functionGP
	    CurrentNetputGP => CurrentFunctionGP%FirstNetputGP
	    NumNetputGPDone = 0
	    do while(NumNetputGPDone.lt.CurrentFunctionGP.NumIFields+CurrentFunctionGP.NumOFields)

	      if(CurrentNetputGP%NumTaxTypes.gt.0) then
	        CurrentTaxGP => CurrentNetputGP%FirstTaxGP
	        CurrentTaxMPSGE => CurrentTaxGP%taxlistMPSGE
	        NumTaxDone = 0
	        do while(NumTaxDone.lt.CurrentNetputGP%NumTaxTypes)

		  if (EqualStringsIgnoreCase(ThisConsumer, CurrentTaxMPSGE%a%ident%name)) then
! This tax revenue goes to ThisConsumer. Add tax to revenue
! SUM over relevant sets
	            ThisSubdomain => CurrentTaxGP%TaxRevenueCoeff%arg
	            call WriteSumsGP(ThisSubDomain, SumString, NSum)
	            call CoeffVarRealToGP(CurrentTaxGP%TaxRevenueCoeff, TempString)
! End the sums
	            call EndSumsGP(NSum, EndSumString)
	            if(NumTermsSoFar.gt.0) then
	              StatementStart = '+'
	            else
	              StatementStart = ''
	            end if
	            OutputString = trim(StatementStart)//' '//trim(SumString)// &
		      ' '//trim(TempString)//' '//EndSumString 
! Write this out
	            call WriteOneLineToTAB('    '//OutputString)
! Increment NumTermsSoFar
	    	    NumTermsSoFar = NumTermsSoFar + 1

	          end if

	          NumTaxDone = NumTaxDone + 1
	          CurrentTaxGP => CurrentTaxGP%next
	          CurrentTaxMPSGE => CurrentTaxMPSGE%next
	        end do ! while

!	       End "if(CurrentNetputGP%NumTaxTypes.gt.0) then"
	      end if

	      NumNetputGPDone = NumNetputGPDone + 1
	      CurrentNetputGP => CurrentNetputGP%next
	    end do  ! while(NumNetputGPDone.lt.CurrentFunctionGP.NumIFields+CurrentFunctionGP.NumOFields)

!	   End "if(CurrentFunctionGP%IsProdFunctionGP) then"
	  end if

	  CurrentFunctionGP => CurrentFunctionGP%next
	end do ! CurrentFunction

! End this statement
	call WriteOneLineToTAB('       ;')

	return
	end subroutine DoEquationsOneDemandFunction

	subroutine DoMarketClearingEquations

! Local
	type (endowmentGP), pointer :: CurrentEndowmentGP
	type (domain), pointer :: ThisSubDomain
	type (functionGP), pointer :: CurrentFunctionGP
	type (netputGP), pointer :: CurrentNetputGP
	integer :: kpass, NumTermsSoFar, NSum, NumNetputGPDone, NumEndowmentGPDone
	character(len=512) OutputString, StatementStart, TempString, OutputStringAlls
	character(len=512) SumString, EndSumString, TempString2, TempString3, TempString4
	character(len=512) ConditionalSumString, TempString5
	type (GPCoeffVarReal) :: CurrentSectorPriceVar
	type (symbol), pointer :: CurrentSymbol
	type (field), pointer :: CurrentDemand
	type (domain), pointer :: FoundArgs
	logical :: UseThis, InDField, ProdFound, FoundWalrasEq, Found, ThisIsWalrasEq
	integer :: NumProdFound
	character :: FoundUserSystemType
	character (len=3) :: ch3
! If find commodity in d: field of a $demand function, 
!  record IncomeVar, TotalIncomeCoeff from this $demand function.
! Also record ActivityLevelVar from the $prod: block containing this commodity in a o: line
	type (GPCoeffVarReal) :: IncomeVarDField !! %-change in value of income for this agent
	type (GPCoeffVarReal) :: TotalIncomeCoeffDField !! Total income for this consumer
	type (GPCoeffVarReal) :: ActivityLevelVarDField !! Activity level in the (unique?) $prod block
			! in which this appears in an o: field
	type (GPCoeffVarReal) :: ThisEndowCoeff, ThisEndowVar
	character(len=llst), allocatable, dimension(:) :: AllSets
	character(len=llindex), allocatable, dimension(:) :: AllIndexes
	integer DimAllArgs
	logical :: UseIndexStack
	type(domain), pointer :: CurrentLineDomain, CurrentDeclaredDomain
	integer :: CurrentArgNumber, assno
! For assertion that total value of supply equals total value of demand
	type (GPCoeffVarReal) :: TotalValueCoeff, TotalValueDemandCoeff, TotalValueSupplyCoeff
	character :: DemandSupplyChar, ISChar



! Run through the commodities (the prices). For each one, there is a single
!  linear equation which says demand = supply. 
!  See which $prod and which $demand functions each one is in. 
! Supplies are in o: lines of $prod functions and in e: lines or $demand functions
! Supplies are in i: lines of $prod functions and in d: lines or $demand functions

	write(*,*)'Beginning Market clearing equations'

	FoundWalrasEq = .false.

! [31/03/04] UseIndexStack (or something similar) is necessary in some cases,
!  though not in all. For the present, always do that. Later may test and
!  only do it when required.
!	UseIndexStack = .false. 
	UseIndexStack = .true. 

! Write headings for this section of the TAB file
	call WriteOneLineToTAB(' ')
	OutputString = '![[! Market clearing equations (one for each $commodity)'
	call WriteOneLineToTAB(OutputString)
	OutputString = &
	    '  But, noting Walras law, one of these equations must be omitted !]]!'
	call WriteOneLineToTAB(OutputString)

! Go through the commodities. Find them by going through idlist.

	CurrentSymbol => idlist
	do while(associated(CurrentSymbol))
	  if(CurrentSymbol%type.eq.'commodity') then

	    ThisIsWalrasEq = .false.

! Can we omit this to satisfy Walras law?
! For the present, omit first one which has no arguments (if is one)
	    if(.not.FoundWalrasEq) then
	      call IsNameInGPNameList(CurrentSymbol%name, Found, FoundUserSystemType, FoundArgs)
	      if(Found) then
		if(.not.associated(FoundArgs)) then
! Omit this to satisfy Walras law
		  FoundWalrasEq = .true.
		  ThisIsWalrasEq = .true.

	          call WriteOneLineToTAB(' ')
		  OutputString = '! [Omitting market clearing equation for '//trim(CurrentSymbol%name)
	          call WriteOneLineToTAB(OutputString)
		  OutputString = '   to take account of Walras law.] !'
	          call WriteOneLineToTAB(OutputString)
		end if
	      else
	      end if
	    end if

	    if(ThisIsWalrasEq) then
! Declare WalrasSlackVar
	      call DeclareGPCoeffVar(WalrasSlackVar)
	    end if

! Write Market clearing equation for this commodity.

! First a comment
	    if(.not.ThisIsWalrasEq) then
	      call WriteOneLineToTAB(' ')
	      OutputString = '! Market clearing for '//trim(CurrentSymbol%name)//' !'
	      call WriteOneLineToTAB(OutputString)
	    end if

	    if(MixedTABFile) then
	      StatementStart = 'Equation (Linear)'
	    else
	      StatementStart = 'Equation'
	    end if
! Write the equation differently if this is WalrasEq
	    if(.not.ThisIsWalrasEq) then
	      OutputString = trim(StatementStart)//' E_'//CurrentSymbol%name
	    else if(ThisIsWalrasEq) then
	      OutputString = trim(StatementStart)//' E_'//WalrasSlackVar%name
	    end if
! Write this out
	    write(*,*)'Beginning: '//trim(OutputString)
	    call WriteOneLineToTAB(OutputString)
! Write label
	    if(.not.ThisIsWalrasEq) then
	      call PrepareLabelInfoForTABFile( &
	        'Market clearing equation for commodity '//trim(CurrentSymbol%name), &
		  OutputString)
	    else if(ThisIsWalrasEq) then
	      call PrepareLabelInfoForTABFile(WalrasSlackVar%labelinfo, OutputString)
	    end if
	    call WriteOneLineToTAB(' '//OutputString)

	    if(ThisIsWalrasEq) then

! WalrasEq case. If the omitted equation is "a+b=c+d", we are writing it as
!    "walrasslack = a+b-c-d"
! Write the LHS and "=" here

	      call CoeffVarRealToGP(WalrasSlackVar, TempString)
	      call WriteOneLineToTAB(' '//trim(TempString)//' =')

	    end if ! if(ThisIsWalrasEq) then

! Quantifiers - always range over sets this price was declared with in $commodity section
! [Checked much earlier in the MPSGE module that the price at the start of the relevant
!   d:, e:, i: or o: field was declared appropriately in $commodity section.]
! [31/03/04] Need to get special all indexes since the expressions can be
!   complicated when sets over which price ranges in field is different from
!   sets over which it is declared in the $commodities part
!
	    if(.not.UseIndexStack) then  
	      call ArgToAllsGP(CurrentSymbol%arg, OutputStringAlls)
	    else if(UseIndexStack) then
	      DimAllArgs = GetArgDim(CurrentSymbol%arg)
! If these arrays are allocated, first deallocate them
	      if(allocated(AllSets)) then
	        deallocate(AllSets, AllIndexes)
	      end if
!  Allocate AllSets and AllIndexes arrays
	      allocate(AllSets(DimAllArgs), AllIndexes(DimAllArgs))
	      call ArgToAllsGPUseIndexStack(CurrentSymbol%arg, DimAllArgs, &
		OutputStringAlls, AllSets, AllIndexes)
	    end if


! Marking this endogenous is done in two possible cases below. [At this stage don't know
!   a GPCoeffVarReal whose name is same as CurrentSymbol%name, which prevents us
!    calling routine MarkAsEndogenous now.]

! Go through the $demand functions to see if this commodity appears in
!  a d: field 

	    InDField = .false.
	    CurrentFunctionGP => FunctionGPList

! [No need to continue if have found it in a d: field in a $demand]
	    do while(associated(CurrentFunctionGP).and.(.not.InDField))

	      if (CurrentFunctionGP%IsDemandFunctionGP) then

! Look at the d: field (if is one)
		CurrentDemand => CurrentFunctionGP%functionMPSGE%demand
		if(associated(CurrentDemand)) then
		  InDField = EqualStringsIgnoreCase( &
		              CurrentDemand%identtext, CurrentSymbol%name)
		  if(InDField) then
		    IncomeVarDField = CurrentFunctionGP%IncomeVar
		    TotalIncomeCoeffDField = CurrentFunctionGP%TotalIncomeCoeff
		  end if 

!	         End "if(associated(CurrentDemand)) then"
		end if

!		 End "if (CurrentFunctionGP%IsDemandFunctionGP) then"
	      end if

	      CurrentFunctionGP => CurrentFunctionGP%next

	    end do ! while(associated(CurrentFunctionGP).and.(.not.InDField))

	   if(InDField) then

! The equation says IncomeVarDField = CurrentSymbol%name + ActivityLevelVar
!  where ActivityLevelVar is the activity level in the (unique?) $prod
!   function containing a line o:CurrentSymbol%name
! Go through the $prod blocks to find the relevant one 

	      write(*,*)'Beginning InDField case'
	      ProdFound = .false.
	      NumProdFound = 0
	      CurrentFunctionGP => FunctionGPList
	      do while(associated(CurrentFunctionGP))

		if (CurrentFunctionGP%IsProdFunctionGP) then

! go through its o: fields
		  NumNetputGPDone = 0
		  CurrentNetputGP => CurrentFunctionGP%FirstNetputGP

! [Probably no need to continue if have found it in an o: field in a $prod. But here
!    continue in case this happens in more than one $prod, which is internal error at present.]

		  do while(NumNetputGPDone.lt.CurrentFunctionGP%NumIFields+CurrentFunctionGP%NumOFields)

! If this is an o: field, look at its commodity
		    if(CurrentNetputGP%netputMPSGE%output) then
		      ProdFound = EqualStringsIgnoreCase( &
		              CurrentNetputGP%IOPriceVar%name, CurrentSymbol%name)
		      if(ProdFound) then
			NumProdFound = NumProdFound + 1
			ActivityLevelVarDField = CurrentFunctionGP%ActivityLevelVar
! CurrentSymbol%name must be the name of CurrentNetputGP%IOPriceVar. Record this for writing the equation
			CurrentSectorPriceVar = CurrentNetputGP%IOPriceVar
		      end if

		    end if ! if(CurrentNetputGP%output) then

! Increment NumNetputGPDOne
		    NumNetputGPDone = NumNetputGPDone + 1

	            CurrentNetputGP => CurrentNetputGP%next

		  end do ! while(NumNetputGPDone.lt.CurrentFunctionGP%NumIFields+CurrentFunctionGP%NumOFields)

!		 End "if (CurrentFunctionGP%IsProdFunctionGP) then"
		end if

		CurrentFunctionGP => CurrentFunctionGP%next
	      end do ! while(associated(CurrentFunctionGP))

	      if(NumProdFound.eq.1) then

! Mark this variable as endogenous
	        if(.not.ThisIsWalrasEq) then
		  write(*,*)'Marking as endogenous [InDField case]:', trim(CurrentSectorPriceVar%name)
	          call MarkAsEndogenous(CurrentSectorPriceVar)
		else if(ThisIsWalrasEq) then
		  write(*,*)'Marking as endogenous [InDField case]:', trim(WalrasSlackVar%name)
	          call MarkAsEndogenous(WalrasSlackVar)
		end if
! The equation says IncomeVarDField = ActivityLevelVarDField + CurrentSectorPriceVar
! [Assume that all these have the same arguments - should check this.]
! Quantifiers have been done earlier.
! Never need to use index stack here, I think
!.NO	        if(.not.UseIndexStack) then
                  call CoeffVarRealToGP(IncomeVarDField, TempString)
                  call CoeffVarRealToGP(ActivityLevelVarDField, TempString2)
                  call CoeffVarRealToGP(CurrentSectorPriceVar, TempString3)
!.NO		else if(UseIndexStack) then
!.NO                  call CoeffVarRealToGPUseIndexStack(IncomeVarDField, DimAllArgs, AllSets, AllIndexes, TempString)
!.NO                  call CoeffVarRealToGPUseIndexStack(ActivityLevelVarDField, DimAllArgs, AllSets, AllIndexes, TempString2)
!.NO                  call CoeffVarRealToGPUseIndexStack(CurrentSectorPriceVar, DimAllArgs, AllSets, AllIndexes, TempString3)
!.NO	        end if
! Write the equation
	        if(.not.ThisIsWalrasEq) then
		  OutputString = trim(OutputStringAlls)//' '//trim(TempString)//' = '// &
		    trim(TempString2)//' + '//trim(TempString3)//' ;'
		else if(ThisIsWalrasEq) then
! LHS (walrasslack)and "=" have already been written. There should be no quantifiers so omit OutputStringAlls
		  OutputString = trim(TempString2)//' + '//trim(TempString3)//' - '//trim(TempString)//' ;'
		end if
	        call WriteOneLineToTAB('    '//OutputString)

	      else if(NumProdFound.ne.1) then
! This is an internal error
		write(ch3, '(I3)') NumProdFound
		call InternalErrorStop( &
		  'Found commodity '//trim(CurrentSymbol%name)//' in '//ch3//' o: fields'// &
		  '\n  but expected to find it in exactly one o: field.'// &
		  '\n   [This commodity appears in d: field in a $demand block.]' )
	      end if
	      write(*,*)'Ended InDField case'

!	    End "if(InDField) then"
	   else

	      write(*,*)'Beginning .not.InDField case'
! go through the $prod and $demand functions.
! Do demands when KPASS=1, supply when KPASS=2.
! Do Total Value of demand when KPASS=3
! Do Total Value of supply when KPASS=4
	

	    do kpass=1,4

! When KPASS=3,4 make up name of Coefficient to hold total value of demand or supply
	      if(kpass.eq.3.or.kpass.eq.4) then

! First a comment
	        if(kpass.eq.3) then
	          call WriteOneLineToTAB(' ')
	          OutputString = '! Calculation of Total Value of Demand and Supply for '//trim(CurrentSymbol%name)//' !'
	          call WriteOneLineToTAB(OutputString)
		end if

	        call InitializeGPCoeffVarReal('c', TotalValueCoeff)
		if(kpass.eq.3) then
		  DemandSupplyChar = 'd'
		else
		  DemandSupplyChar = 's'
		end if
	        TotalValueCoeff%name = ConstructGPName('c', 's', CurrentSymbol%arg, &
					'TV'//DemandSupplyChar//'_'//trim(CurrentSymbol%name), 'TV'//DemandSupplyChar)
	        call upper(TotalValueCoeff%name)   ! Convert name to upper case
	        TotalValueCoeff%arg => CurrentSymbol%arg
		if(kpass.eq.3) then 
	          TotalValueCoeff%labelinfo = 'Total value of demand for commodity '//CurrentSymbol%name
		else if(kpass.eq.4) then
	          TotalValueCoeff%labelinfo = 'Total value of supply for commodity '//CurrentSymbol%name
		end if
		if(kpass.eq.3) then
		  TotalValueDemandCoeff = TotalValueCoeff
		else if(kpass.eq.4) then
		  TotalValueSupplyCoeff = TotalValueCoeff
		end if
! Declare this
	        call DeclareGPCoeffVar(TotalValueCoeff)

! Start a formula(always) for this coeff
		StatementStart = 'Formula (Always)'
! Write this out
	        write(*,*)'Beginning: '//trim(StatementStart)
	        call WriteOneLineToTAB(StatementStart)
! Labelling info for this formula
	        call PrepareLabelInfoForTABFile( TotalValueCoeff%labelinfo, OutputString)
	        call WriteOneLineToTAB(' '//OutputString)
! quantifiers and LHS and "="
		if(.not.UseIndexStack) then
	          call ArgToAllsGP(CurrentSymbol%arg, OutputStringAlls)
	          call CoeffVarRealToGP(TotalValueCoeff, TempString)
		else if(UseIndexStack) then
	          call ArgToAllsGPUseIndexStack(CurrentSymbol%arg, DimAllArgs, &
		    OutputStringAlls, AllSets, AllIndexes)
	          call CoeffVarRealToGPUseIndexStack( &
				TotalValueCoeff, &
			        DimAllArgs, AllSets, AllIndexes, TempString)
		end if
! Write this 
		OutputString = trim(OutputStringAlls)//' '//trim(TempString)//' ='
	        call WriteOneLineToTAB('   '//OutputString)

	      end if ! if(kpass.eq.3.or.kpass.eq.4) then

! NumTermsSoFar is used to decide if need "+" before this term
	      NumTermsSoFar=0

	      CurrentFunctionGP => FunctionGPList
	      do while(associated(CurrentFunctionGP))

		if (CurrentFunctionGP%IsDemandFunctionGP) then

! d: field has been checked earlier in this routine. So nothing to do here if KPASS=1.


		  if(kpass.eq.2.or.kpass.eq.4) then

! go through the e: fields if this is a $demand and KPASS=2
		    CurrentEndowmentGP => CurrentFunctionGP%FirstEndowmentGP
		    NumEndowmentGPDone = 0
		    do while(NumEndowmentGPDone.lt.CurrentFunctionGP.NumEFields)
		      UseThis = EqualStringsIgnoreCase( &
		          CurrentEndowmentGP%endowmentMPSGE%f%identtext, &
			  CurrentSymbol%name)

		      if(UseThis) then

! Sum over all sets in EndowQFieldValueCoeff which are not arguments of this commodity
! Use subroutine findsubdomain to work this out.
!       subroutine findsubdomain(d1,d2,subdomain) determines those elements of d1 which are not in d2.

		        call findsubdomain(CurrentEndowmentGP%EndowQFieldValueCoeff%arg, &
			      CurrentEndowmentGP%EndowPriceVar%arg,ThisSubdomain)
	                call WriteSumsGP(ThisSubDomain, SumString, NSum)
! If set price ranges over in e: line (or o: line or i: line) is different from the
!  set that argument of the price ranges over (as declared in $commodities section),
!   need to add a conditional sum.
			CurrentLineDomain => CurrentEndowmentGP%EndowPriceVar%arg
			CurrentDeclaredDomain => CurrentSymbol%arg
			CurrentArgNumber = 1
			do while(associated(CurrentLineDomain))
			  if(.not.EqualStringsIgnoreCase(CurrentLineDomain%s%name, CurrentDeclaredDomain%s%name)) then
! [Even do this when set is an element]
!.NO	                    if(INDEX(CurrentLineDomain%s%name,'"').eq.0) then
	                      call WriteConditionalSumGP(CurrentLineDomain, &
				CurrentArgNumber, DimAllArgs, AllSets, AllIndexes, ConditionalSumString)
! Add this to SumString and increase NSum
			      SumString = trim(SumString)//ConditionalSumString
			      NSum = NSum + 1
!.NO			    end if
			  end if ! if(.not.EqualStringsIgnoreCase(CurrentDomain%s%name, DeclaredDomain%s%name)) then
			  CurrentLineDomain => CurrentLineDomain%next
			  CurrentDeclaredDomain => CurrentDeclaredDomain%next
			  CurrentArgNumber = CurrentArgNumber + 1
			end do ! while(associated(CurrentLineDomain))

!	  type (GPCoeffVarReal) :: EndowQFieldQuantityVar  !! %-change in quantity of each endowment  
!	  type (GPCoeffVarReal) :: EndowQFieldValueCoeff  !! value (level) of endowment
! SUM  of terms: EndowQFieldValueCoeff * EndowQFieldQuantityVar

			ThisEndowCoeff = CurrentEndowmentGP%EndowQFieldValueCoeff
			ThisEndowVar = CurrentEndowmentGP%EndowQFieldQuantityVar 

			if(.not.UseIndexStack) then
	                  call CoeffVarRealToGP(ThisEndowCoeff, TempString)
	                  call CoeffVarRealToGP(ThisEndowVar, TempString2)
			else if(UseIndexStack) then
	                  call CoeffVarRealToGPUseIndexStack( &
				ThisEndowCoeff, &
			        DimAllArgs, AllSets, AllIndexes, TempString)
	                  call CoeffVarRealToGPUseIndexStack( &
				ThisEndowVar, &
				DimAllArgs, AllSets, AllIndexes, TempString2)
			end if

! End the sums
	                call EndSumsGP(NSum, EndSumString)
			if(ThisIsWalrasEq.and.kpass.eq.2) then
! WalrasEq case. If the omitted equation is "a+b=c+d", we are writing it as
!    "walrasslack = a+b-c-d" 
!  Hence need "-" here when kpass=2
		          OutputString='-'
		        else if(NumTermsSoFar.gt.0) then
		          OutputString='+'
		        else
		          OutputString=''
		        end if
			if(kpass.eq.2) then
! The equation case
		          OutputString = trim(OutputString)//' '//trim(SumString)//' '// &
			    trim(TempString)//' * '//trim(TempString2)//' '//trim(EndSumString)
			else if(kpass.eq.4) then
! The formula case
		          OutputString = trim(OutputString)//' '//trim(SumString)//' '// &
			    trim(TempString)//' '//trim(EndSumString)
			end if
! Write this line
	    	        call WriteOneLineToTAB('     '//OutputString)

! Increment NumTermsSoFar
		        NumTermsSoFar = NumTermsSoFar + 1

		      end if  ! if(UseThis) then

		      NumEndowmentGPDone = NumEndowmentGPDone + 1
		      CurrentEndowmentGP => CurrentEndowmentGP%next
		    end do  ! while(NumEndowmentGPDone.lt.CurrentFunctionGP.NumEFields)

!		   End "else if(kpass.eq.2.or.kpass.eq.4) then"
		  end if


!		 End "if (CurrentFunctionGP%IsDemandFunctionGP) then"
		else if (CurrentFunctionGP%IsProdFunctionGP) then

! go through the relevant netputGPs if this is a $prod function

		  CurrentNetputGP => CurrentFunctionGP%FirstNetputGP
		  NumNetputGPDone = 0
		  do while(NumNetputGPDOne.lt.(CurrentFunctionGP%NumIFields+CurrentFunctionGP%NumOFields))
		    UseThis=((kpass.eq.2.or.kpass.eq.4).and.CurrentNetputGP%netputMPSGE%output).or. &
		            ((kpass.eq.1.or.kpass.eq.3).and.(.not.CurrentNetputGP%netputMPSGE%output))
		    UseThis = UseThis .and. & 
		      EqualStringsIgnoreCase( &
		        CurrentNetputGP%netputMPSGE%f%identtext, &
			CurrentSymbol%name)

		    if(UseThis) then

! Output Quantifiers if KPASS=1 and NumTermsSoFar=0. [Have to wait until now since
!		 cannot tell the args of this from info about is as a symbol. So wait
!		  until find it as the price at start of o: or d: line.]

! Tricky example from SJ.MPS
! $commodities:
!  pc(SECT)              ! price of commodity
! $prod:xcom(SECT)          
!* Tricky. Can infer that SECT1 is a subset of SECT here since pc is defined over SECT
!*  but is also used over SECT1.
!  o:pc(SECT)         q:DVCOM(SECT)
!  i:pc(SECT1)        q:DVCOMIN(SECT1,SECT)
!  i:pf(FAC)          q:DVFACIN(FAC,SECT)
! Here the ALLs in market clearing equation for pc must be over the set over which pc 
!  was declared in the $commodities section
! And must map indices accordingly. That is, need (is this correct??)
! (all,sect_1,SECT) DVCOM(sect_1)*pc(sect_1) = SUM(sect1_1,SECT1, DVCOMIN(sect1_1,sect_1)*pc(sect1_1)
! Considerable care will be needed to map the indices correctly. **

		      if(kpass.eq.1.and.NumTermsSoFar.eq.0) then
! [Quantifiers are done earlier]
! Write this line
	                call WriteOneLineToTAB(' '//OutputStringAlls)
! Mark this variable as endogenous
	                call MarkAsEndogenous(CurrentNetputGP%IOPriceVar)
		      end if

! Sum over all sets in IOQFieldValueCoeff which are not arguments of this commodity
! Use subroutine findsubdomain to work this out.
!       subroutine findsubdomain(d1,d2,subdomain) determines those elements of d1 which are not in d2.

		      call findsubdomain(CurrentNetputGP%IOQFieldValueCoeff%arg, &
			CurrentNetputGP%IOPriceVar%arg,ThisSubdomain)
	              call WriteSumsGP(ThisSubDomain, SumString, NSum)

! If set price ranges over in e: line (or o: line or i: line) is different from the
!  set that argument of the price ranges over (as declared in $commodities section),
!   need to add a conditional sum.
		      CurrentLineDomain => CurrentNetputGP%IOPriceVar%arg
		      CurrentDeclaredDomain => CurrentSymbol%arg
		      CurrentArgNumber = 1
		      do while(associated(CurrentLineDomain))
			  if(.not.EqualStringsIgnoreCase(CurrentLineDomain%s%name, CurrentDeclaredDomain%s%name)) then
! [Even do this when set is an element]
!.NO	                    if(INDEX(CurrentLineDomain%s%name,'"').eq.0) then
	                      call WriteConditionalSumGP(CurrentLineDomain, &
				CurrentArgNumber, DimAllArgs, AllSets, AllIndexes, ConditionalSumString)
! Add this to SumString and increase NSum
			      SumString = trim(SumString)//ConditionalSumString
			      NSum = NSum + 1
!.NO			    end if
			  end if ! if(.not.EqualStringsIgnoreCase(CurrentDomain%s%name, DeclaredDomain%s%name)) then
			  CurrentLineDomain => CurrentLineDomain%next
			  CurrentDeclaredDomain => CurrentDeclaredDomain%next
			  CurrentArgNumber = CurrentArgNumber + 1
		      end do ! while(associated(CurrentLineDomain))

! SUM  of terms: IOQFieldValueCoeff * IOQFieldQuantityVar
		      if(.not.UseIndexStack) then
	                call CoeffVarRealToGP(CurrentNetputGP%IOQFieldValueCoeff, TempString)
	                call CoeffVarRealToGP(CurrentNetputGP%IOQFieldQuantityVar, TempString2)
		      else if(UseIndexStack) then
	                call CoeffVarRealToGPUseIndexStack( &
			    CurrentNetputGP%IOQFieldValueCoeff, &
			    DimAllArgs, AllSets, AllIndexes, TempString)
	                call CoeffVarRealToGPUseIndexStack( &
			    CurrentNetputGP%IOQFieldQuantityVar, &
			    DimAllArgs, AllSets, AllIndexes, TempString2)
		      end if
! End the sums
	              call EndSumsGP(NSum, EndSumString)
		      if(NumTermsSoFar.gt.0) then
		        OutputString='+'
		      else
		        OutputString=''
		      end if
		      if(kpass.le.2) then
! The equation case
		        OutputString = trim(OutputString)//' '//trim(SumString)//' '// &
			  trim(TempString)//' * '//trim(TempString2)//' '//trim(EndSumString)
		      else if(kpass.ge.3) then
! The formula case
		        OutputString = trim(OutputString)//' '//trim(SumString)//' '// &
			  trim(TempString)//' '//trim(EndSumString)
		      end if
! Write this line
		      write(*,*)'Writing line: '//trim(OutputString),' to TAB file'
	    	      call WriteOneLineToTAB('     '//OutputString)

! Increment NumTermsSoFar
		      NumTermsSoFar = NumTermsSoFar + 1

!		     End "if(UseThis) then"
		    end if

! Increment NumNetputGPDOne
		    NumNetputGPDone = NumNetputGPDone + 1

		    CurrentNetputGP => CurrentNetputGP%next
		  end do ! while(associated(CurrentNetputGP))


!		 End "else if (CurrentFunctionGP%IsProdFunctionGP) then"
		end if

	        CurrentFunctionGP => CurrentFunctionGP%next
	      end do !  while(associated(CurrentFunctionGP))

	      if(kpass.eq.1) then
! "=" sign between the two sides
! WalrasEq case. If the omitted equation is "a+b=c+d", we are writing it as
!    "walrasslack = a+b-c-d"
! Hence nothing here if this is WalrasEq.
	        if(.not.ThisIsWalrasEq) then 
		  call WriteOneLineToTAB('             =')
		end if
	      else if(kpass.eq.2) then
! ";" to end equation
		call WriteOneLineToTAB('             ;')
	      else if(kpass.eq.3.or.kpass.eq.4) then
! ";" to end formula
		call WriteOneLineToTAB('             ;')
	      end if

	      if(kpass.eq.4) then

! Add assertions about value of demand and supply
! Do slightly stricter assertion for intial data, less strict at other steps

! First a comment
	        call WriteOneLineToTAB(' ')
	        OutputString = '! Checks that Total Value of Supply and Demand are equal for '//trim(CurrentSymbol%name)//' !'
	        call WriteOneLineToTAB(OutputString)

	        do assno=1,2

! Do Assertion (Inital) when KPASS=1	
! Do Assertion (Always) when KPASS=2

	          if(assno.eq.1) then
	            ISChar = 'I'
	          else if(assno.eq.2) then
	            ISChar = 'S'
	          end if	

! Assertion that total value of demand = total value of supply
	          if(assno.eq.1) then
	    	    StatementStart = 'Assertion (Initial)'
	  	  else if(assno.eq.2) then
	    	    StatementStart = 'Assertion (Always)'
	  	  end if
	          OutputString = trim(StatementStart)//&
		    ' # Total Value of Demand = Total Value of Supply for '//&
		    trim(CurrentSymbol%name)//' #'
! Write this out
	  	  call WriteOneLineToTAB(OutputString)
! Quantifiers
	  	  call ArgToAllsGP(CurrentSymbol%arg, OutputString)
! Write this out
	  	  call WriteOneLineToTAB(OutputString)
! LHS
	          call CoeffVarRealToGP(TotalValueDemandCoeff, TempString)
! RHS
	          call CoeffVarRealToGP(TotalValueSupplyCoeff, TempString3)

! Cannot expect exact equality. So instead write it as
!  Demand > 0.9999*Supply and Demand < 1.0001*Supply
! [Equality]	OutputString = trim(OutputString)//' '//trim(TempString)// &
! [Equality]	  ' = '//trim(TempString3)//' ;'
! Use BALTOLI or BALTOLS
!	  OutputString = trim(OutputString)//' '//trim(TempString)// &
!	    ' GT 0.9999*'//trim(TempString3)//' AND '//trim(TempString)// &
!	    ' LT 1.0001*'//trim(TempString3)//' ;'

	          TempString4 = '['//trim(TempString)// ' GT (1-BALTOL'//ISChar//')*'//trim(TempString3)//'] AND'
	          write(*,*)'TempString4 = ', trim(TempString4)
	          TempString2 = '   '//trim(TempString4)
	  	  call WriteOneLineToTAB(TempString2)
	  	  TempString5 = '['//trim(TempString)// ' LT (1+BALTOL'//ISChar//')*'//trim(TempString3)//'] ;'
	  	  write(*,*)'TempString5 = ', trim(TempString5)
	  	  TempString2 = '   '//trim(TempString5)
	  	  call WriteOneLineToTAB(TempString2)
		end do ! do assno=1,2

	      end if ! if(kpass.eq.4) then


!	     End "do kpass=1,4"
	    end do

	      write(*,*)'Ended .not.InDField case'

!	    End "if(InDField) then..else.."
	   end if

!	   End "if(CurrentSymbol%type.eq.'commodity') then"
	  end if

	  CurrentSymbol => CurrentSymbol%next
	end do 	! while(associated(CurrentSymbol))

	call WriteOneLineToTAB(' ')
	OutputString = '! End of market clearing equations !'
	call WriteOneLineToTAB(OutputString)

	write(*,*)'Completed Market clearing equations'


	return
	end subroutine DoMarketClearingEquations

	subroutine StartTABGP

! Local
	type (symbol), pointer :: ThisSymbol
	character (len=llhd) HeadString
	character (len=80) TextLine, OutputString
	character (len=llst) :: BigSetName, SmallSetName, SetName, SetName2, SmallSetName2
	integer :: NumSubsetsHere, s1, s2, n1

10005	format(a)

! Coefficients used for checking balance
	call WriteOneLineToTAB(' ')
	OutputString = '! Coefficients used to check balance of data !'
	call WriteOneLineToTAB(OutputString)
	OutputString = 'Coefficient (Parameter) BALTOLI'
	call WriteOneLineToTAB(OutputString)
	OutputString = '  # Used to check initial database is balanced # ;'
	call WriteOneLineToTAB(OutputString)
	OutputString = 'Coefficient (Parameter) BALTOLS'
	call WriteOneLineToTAB(OutputString)
	OutputString = '  # Used to check database stays balanced at all steps # ;'
	call WriteOneLineToTAB(OutputString)
	OutputString = 'Formula (Initial) BALTOLI = 0.0001 ;' 
	call WriteOneLineToTAB(OutputString)
	OutputString = 'Formula (Initial) BALTOLS = 0.001 ;' 
	call WriteOneLineToTAB(OutputString)

	call WriteOneLineToTAB(' ')
	call WriteOneLineToTAB('File Input # contains data and set elements # ;')

! Read all sets

	ThisSymbol => idlist
	do while(associated(ThisSymbol))
	  if (ThisSymbol%type.eq.'set') then
! Declare this set
	    call GetReadHeader(ThisSymbol%name, HeadString)
	    SetName = trim(ThisSymbol%name)
! In the idlist are sets corresponding to quoted elements (eg "lab") used. Do not
!  want to declare these or read their elements
	    if(INDEX(SetName,'"').eq.0) then
	      call upper(SetName)
	      TextLine='Set '//trim(SetName)//' Read Elements from File Input Header "'//HeadString//'" ;'
	    else
! Here make up a set to hold this single element
	      SetName2 = GetSingletonSetName(SetName)
	      n1=LEN_TRIM(SetName)
	      TextLine='Set '//trim(SetName2)//' ('//SetName(2:n1-1)//') ;'
	    end if
	    call WriteOneLineToTAB(TextLine)
	  end if
	  ThisSymbol => ThisSymbol%next
	end do

! Write the subset statements (if any)
! First see if are any - if so, write heading

	NumSubsetsHere = 0
	do s1=1,NumSets
	  do s2=1,NumSets
	    if(IsSubsetGP(s1,s2)) then
	        NumSubsetsHere = NumSubsetsHere + 1
	    end if
	  end do
	end do

	if(NumSubsetsHere.gt.0) then
	  call WriteOneLineToTAB(' ')
	  TextLine = '! Subset declarations (inferred from MPS file) ! '
	  call WriteOneLineToTAB(TextLine)
	end if

	do s1=1,NumSets
	  do s2=1,NumSets
	    if(IsSubsetGP(s1,s2)) then
! set number s1 is a subset of set number s2
! Don't do this when small set is an element
	      SmallSetName = GetSetName(s1)
	      if(INDEX(GetSetName(s1),'"').eq.0) then
	        SmallSetName2 = GetSetName(s1)
	      else
	        SmallSetName2 = GetSingletonSetName(SmallSetName)
	      end if
	      BigSetName = GetSetName(s2)
	      call upper(BigSetName)
	      TextLine='Subset '//trim(SmallSetName2)//' is subset of '//trim(BigSetName)//' ;'
	      call WriteOneLineToTAB(TextLine)
	    end if
	  end do
	end do

	return
	end subroutine StartTABGP

	subroutine WriteSumsGP(ThisDomain, OutputString, NSum)
! Input
	type (domain), pointer :: ThisDomain
! Output
	character (len=*), intent(out) :: OutputString
	integer, intent(out) :: NSum
! Local
	type (domain), pointer :: CurrentDomain
	character (len=128) :: CurrentIndex, TempString 

! The RHS  "SUM{<index>,<set>,<InputQuantity>} ;"
	NSum = 0
	CurrentDomain => ThisDomain
	OutputString = ''
	do while(associated(CurrentDomain))
! Sum over this set
	  call GPSet2Index(CurrentDomain%s, CurrentIndex)
	  TempString = 'SUM{'//trim(CurrentIndex)//','//trim(CurrentDomain%s%name)//','
	  OutputString = trim(OutputString)//trim(TempString)
	  NSum = NSum + 1
	  CurrentDomain => CurrentDomain%next
	end do

	return
	end subroutine WriteSumsGP

	subroutine WriteConditionalSumGP( &
! Input
	  ThisDomain, CurrentArgNumber, &
	  DimTheseArgs, AllSets, AllIndexes, &
! Output
	  OutputString)
! Input
	type (domain), pointer :: ThisDomain
	integer, intent(in) :: DimTheseArgs, CurrentArgNumber
	character (len=*), intent(in) :: AllSets(DimTheseArgs)
	character (len=*), intent(in) :: AllIndexes(DimTheseArgs)
! Output
	character (len=*), intent(out) :: OutputString
! Local
	type (domain), pointer :: CurrentDomain
	type (symbol), pointer :: CurrentSet
	character (len=128) :: CurrentIndex, TempString
	character (len=llst) :: SumSetName
	integer :: ThisArgNumber 

! "SUM{<index>,<set> : <condition>,"
	OutputString = ''
! Run through ThisDomain to find the set symbol for CurrentArgNumber
	CurrentDomain => ThisDomain
	ThisArgNumber = 0
	do while(ThisArgNumber.lt.CurrentArgNumber)
	  ThisArgNumber = ThisArgNumber + 1
	  if(ThisArgNumber.eq.CurrentArgNumber) then
! Sum over CurrentSet - taken from ThisDomain
	    CurrentSet => CurrentDomain%s	    
	  end if
	  CurrentDomain => CurrentDomain%next
	end do ! while(ThisArgNumber.lt.CurrentArgNumber)
! If set is a quoted element, convert to the corresponding GEMPACK set
	if(INDEX(CurrentSet%name,'"').gt.0) then
	  SumSetName = GetSingletonSetName(CurrentSet%name)
	else
	  SumSetName = CurrentSet%name
	end if
! Sum over this set
	call GPSetName2Index(SumSetName, CurrentIndex)
	TempString = 'SUM{'//trim(CurrentIndex)//','//trim(SumSetName)
	OutputString = trim(OutputString)//trim(TempString)
! Add the condition ": CurrentIndex = <index-from-stack>,"
	TempString = ': '//trim(CurrentIndex)//' = '//trim(AllIndexes(CurrentArgNumber))
	OutputString = trim(OutputString)//trim(TempString)//','
!	write(*,*)'WriteConditionalSumGP. OutputString is: ', trim(OutputString)

	return
	end subroutine WriteConditionalSumGP


	subroutine EndSumsGP(NSum, OutputString)
! Input
	integer, intent(in) :: NSum
! Output
	character (len=*), intent(out) :: OutputString
! Local
	integer :: i
 
	OutputString=''
	do i = 1, NSum
	  OutputString = trim(OutputString)//'}'
	end do

	return
	end subroutine EndSumsGP

	subroutine OutputMPSFunctionAsGPComment(ThisFunction)
! Input
	type (function), pointer :: ThisFunction
! Local
	type (GPCoeffVarReal) :: ThisTopVar
	type (GPCoeffVarReal) :: ThisElasticity
	character (len=128) :: TempString
	type (record), pointer :: CurrentRecord
	type (line), pointer :: CurrentLine, ThisLine
!	integer NRecords, KRecord
	logical :: FoundNextFunction
	integer :: NumLineWritten

10015	format(a)

	TempString = trim(ThisFunction%firstrecord%firstfield%data)
! Blank line
 	write(iotabgp,10015) ' '
 	write(iotabgp,10015) '! '//trim(ThisFunction%firstrecord%firstfield%label)// &
	  ':'//trim(ThisFunction%firstrecord%firstfield%data)//' !'
!	NRecords = ThisFunction%nrecord
!	write(*,*)'NRecords = ', NRecords
!	KRecord = 0
! Write out the records as a strong comment
 	write(iotabgp,10015) '![[!'

	NumLineWritten = 0

! ThisLine may consist of several actual lines from the file
	    ThisLine => ThisFunction%firstrecord%line
	    FoundNextFunction = .false.
	    do while(associated(ThisLine).and. (.not.FoundNextFunction))

! Don't know where this function ends. So test every line to see if it starts
!  with $. If so, set FoundNextFunction true, which will stop the writing.
! But don't test the first line, which is the start of the current function
	      if(NumLineWritten.gt.0.and.ThisLine%text.ne.' ') then
		TempString = ThisLine%text
! Remove blanks at the start until find first non-blank character
		do while(TempString(1:1).eq.' ')
		  TempString = TempString(2:)						      
		end do ! while(TempString(1:1).eq.' ')
		if(TempString(1:1).eq.'$') then
		  FoundNextFunction = .true.
		end if
	      end if ! if(NumLineWritten.gt.0.and.ThisLine%text.ne.' ') then
	      if(.not.FoundNextFunction) then
! Lines may be too long. Use routine OutputOneMPSLineToGP to break them up
! (was)	      write(iotabgp,10015) trim(CurrentLine%text)
!	        write(*,*)'Calling OutputOneMPSLineToGP. ThisLine%text is: ',trim(ThisLine%text) 
	        call OutputOneMPSLineToGP(ThisLine%text)
!	        write(*,*)'Called OutputOneMPSLineToGP.'
		NumLineWritten = NumLineWritten + 1
	      end if ! if(.not.FoundNextFunction) then
	      ThisLine => ThisLine%next
	    end do ! while(associated(ThisLine))

 	write(iotabgp,10015) '!]]!'

	return
	end subroutine OutputMPSFunctionAsGPComment

! Maybe need to make this a recursive subroutine
!??	recursive subroutine OutputOneMPSLineToGP(ThisLine)
	subroutine OutputOneMPSLineToGP(ThisLine)
! Input
	character (len=*), intent(in) :: ThisLine
! Local
	character (len=255) :: TempLine
	integer :: n1
	logical :: Done
!
10015	format(a)
!
	Done = .false.
	if(LEN_TRIM(ThisLine).le.80) then
	  write(iotabgp,10015) trim(ThisLine)
	else
	  TempLine = trim(ThisLine)
! Attempt to shorten the line by replacing two blanks by a single blank.
!   If get short enough, stop
!  [Don't replace two blanks at the start - hence "(2:)" below.]
! ** Must use "trim" here since don't want to be messing with all the trailing blanks
	  do while(.not.Done .and. (INDEX(trim(TempLine(2:)),'  ').gt.0))
	    n1=INDEX(trim(TempLine(2:)), '  ')
! Note that n1 is position in TempLine(2:) so n1+1 is position in TempLine.
	    TempLine = TempLine(1:n1)//' '//TempLine(n1+3:)
	    if(LEN_TRIM(TempLine).le.80) then
	      write(iotabgp,10015) trim(TempLine)
	      Done = .true.
	    end if
	  end do ! while
!	  write(*,*)'OutputOneMPSLineToGP. Finished 2 blanks'
! Here line is still too long. Break it up into chunks
	  do while(.not.Done .and. (LEN_TRIM(TempLine).gt.80))
	    n1=INDEX(TempLine(1:80), ':',back=.true.)
	    if(n1.gt.0) then
!		Output up to the preceding character
	      write(iotabgp,10015) trim(TempLine(1:n1-2))
	      TempLine = TempLine(n1-1:)
	    end if
	    if(LEN_TRIM(TempLine).le.80) then
	      write(iotabgp,10015) trim(TempLine)
	      Done = .true.
	    else
! Here line is longer than 80. Split it arbitrarily
	      n1=INDEX(TempLine(1:80),' ',back=.true.)
	      if(n1.gt.0) then
!		Output up to just before this space
!		write(*,*)'OutputOneMPSLineToGP. Splitting at space in pos ', n1
	        write(iotabgp,10015) trim(TempLine(1:n1-1))
	        TempLine = TempLine(n1:)
	      else
! Desparate here
	        write(iotabgp,10015) trim(TempLine(1:70))
!		write(*,*)'OutputOneMPSLineToGP. Splitting at pos 70'
		TempLine = TempLine(71:)
	      end if
	    end if
	  end do ! while
! Can never get to here??
	  if(.not.Done) then
! Just output what we have. Should not happen I think.
	      write(iotabgp,10015) trim(TempLine)
	  end if	  
	end if

	return
	end subroutine OutputOneMPSLineToGP

	logical function IsReal(ThisString)
! Input
	character (len=*) ThisString
! Local
	integer LengthString, i
	character (len=12) RealChar
	character ThisChar
	intrinsic Index

	RealChar='0123456789.-'

	LengthString = LEN(trim(ThisString))
	IsReal = .true .
	do i = 1, LengthString
	  ThisChar = ThisString(i:i)
	  if(Index(RealChar,ThisChar).eq.0) then
	    IsReal = .false.
	  end if
	end do

	return
	end function IsReal

	subroutine GetProdPriceQuantity(ThisFunction, IOType, ThisPrice, ThisQuantity)
! Input
	type (function), pointer :: ThisFunction
	character, intent(in) :: IOType   ! 'i' or 'o'
! Output
	type (GPCoeffVarReal) :: ThisPrice, ThisQuantity
! Local
	type (record), pointer :: r
	type (field), pointer :: fnf
	type (nest), pointer :: outputs, inputs
	type (netputlist), pointer :: elements
	type (netput), pointer :: ThisNetPut
	type (field), pointer :: IOField, qref, pref
	logical :: UseThis, UseThis1, UseThis2

!	write(*,*)' GetProdPriceQuantity: IOType=', IOType
	r => ThisFunction%firstrecord
	fnf => r%firstfield
	if(IOType.eq.'o') then
	  outputs => ThisFunction%outputs
	  elements => outputs%elements
	else if(IOType.eq.'i') then
	  inputs => ThisFunction%inputs
	  elements => inputs%elements
!	else
!	  error
	end if
	do while (associated(elements)) 
	  ThisNetPut => elements%np
	  IOField => ThisNetPut%f   ! I: or O: field
	  qref => ThisNetPut%qref
	  pref => ThisNetPut%pref
	  UseThis1 = (IOType.eq.'o' .and. ThisNetPut%output)  
	  UseThis2 = (IOType.eq.'i' .and. .not.ThisNetPut%output)
	  UseThis = UseThis1.or.UseThis2  
	  if (UseThis) then
	    ThisPrice%type = 'v' ! Variable
	    ThisPrice%name = IOField%ident%name
	    ThisPrice%arg => IOField%arg	    
	    ThisQuantity%type = 'c' ! Coefficient
	    ThisQuantity%name = qref%ident%name
	    ThisQuantity%arg => qref%arg	    
!  Are done - return
	    return
	  end if
	  elements => elements%next
	end do

	return
	end subroutine GetProdPriceQuantity

	subroutine GetInputsSubDomain(ThisFunction, ThisSubdomain)
! Input
	type (function), pointer :: ThisFunction
! Output
	type (domain), pointer :: ThisSubdomain
! Local
	type (record), pointer :: r
	type (field), pointer :: fnf
	type (nest), pointer :: outputs, inputs
	type (netputlist), pointer :: elements
	type (netput), pointer :: ThisNetput
	logical :: DoneOk

	DoneOk = .false .
	r => ThisFunction%firstrecord
	fnf => r%firstfield
	inputs => ThisFunction%inputs
	elements => inputs%elements
	if (associated(elements)) then 
	  ThisNetput => elements%np
	  if(associated(ThisNetput)) then
	    ThisSubdomain => ThisNetput%sd
	    DoneOk = .true .
	  end if
	end if

	if(.not.DoneOk) then
	  write(*,*)'GetInputsSubDomain - not done ok'
	end if

	return
	end subroutine GetInputsSubDomain

	subroutine ReadCoeffGP(ThisCoefficient)
! Input
	type (GPCoeffVarReal) :: ThisCoefficient
! Local
	character (len=llhd) :: HeadString
	character (len=80) :: TextLine, CoeffName 	

10005	format(a)

	call GetReadHeader(ThisCoefficient%name, HeadString)
	CoeffName=trim(ThisCoefficient%name)
	call upper(CoeffName)
	TextLine = 'Read '//trim(CoeffName)//' from File Input Header "'//HeadString//'" ;'
	write(iotabgp, 10005) trim(TextLine)

	return
	end subroutine ReadCoeffGP	

	subroutine InitializeGPCoeffVarReal(ThisType, ThisGPCoeffVarReal)
! Input
	character (len=*), intent(in) :: ThisType  ! 'v', 'c' or 'r'
! Output
	type (GPCoeffVarReal), intent(out) :: ThisGPCoeffVarReal
!
!  character :: type  ! 'c' for coefficient or 'v' for variable or 'r' for real
!  character (len=15) :: name    ! Coeff/Var name or Real as string
!  type (domain), pointer :: arg      ! Sets over which arguments range
!  character (len=80) :: labelinfo  ! labelling information (between #)
! Next two only apply to a variable
!  character :: vartype  ! 'p' for %-change, 'c' for change
!  type (GPCoeffVarReal) :: OrigLevel  ! Coeff or Real which is the ORIG_LEVEL ! ignore if OrigLevel%name is ' '
! Next only applies to a coefficient
!  logical :: IsParam  ! Tells if this coefficient is a parameter
!
	ThisGPCoeffVarReal%type = ThisType
	ThisGPCoeffVarReal%name = ' '
	nullify(ThisGPCoeffVarReal%arg)
	ThisGPCoeffVarReal%labelinfo = ' '
	if(ThisType.eq.'v') then
	  ThisGPCoeffVarReal%vartype = 'p'
	else
	  ThisGPCoeffVarReal%vartype = ' '
	end if
	ThisGPCoeffVarReal%OrigLevel = ' '
	ThisGPCoeffVarReal%IsParam = .false.

	return
	end subroutine InitializeGPCoeffVarReal

	subroutine MarkAsEndogenous(ThisGPVar)
! Input/Output
	type (GPCoeffVarReal) :: ThisGPVar
! Local
	type (GPName), pointer :: CurrentGPName
	logical :: found

! ?? Later may need to take into account sets??

! Run through GPNameList to find this name, and mark
	found = .false.
	CurrentGPName => GPNameList
	do while(associated(CurrentGPName).and.(.not.found))
	  if(EqualStringsIgnoreCase(ThisGPVar%name,CurrentGPName%name)) then
	    found = .true.
	    CurrentGPName%ExogEndogStatus = 'n'
	  end if
	  CurrentGPName => CurrentGPName%next
	end do ! while(associated(CurrentGPName).and.(.not.found))	

	return
	end subroutine MarkAsEndogenous

	subroutine WriteArguments(FirstArg, OutputText, LenOutput)
! Input
	type (domain), pointer :: FirstArg
! Output
	character (len=255) OutputText
	integer LenOutput
! Local
	type (domain), pointer :: CurrentArg
	type (symbol), pointer :: CurrentSet

10017   format(2x,3a)
	OutputText = ''
	LenOutput = 0	
	CurrentArg => FirstArg	
	do while(associated(CurrentArg))
	  CurrentSet => CurrentArg%s
	  OutputText(LenOutput+1:) = trim(CurrentSet%name)//' '
	  LenOutput = LenOutput + LEN(trim(CurrentSet%name))+1
	  CurrentArg => CurrentArg%next
	end do
!	write(*,10017) 'Arguments: ', OutputText(1:LenOutput) 

	return
	end subroutine WriteArguments 
	
	subroutine WriteRecord(DescText, ThisRecord)
! Input
	character(len=*), intent(in) :: DescText
	type (record), pointer :: ThisRecord
! Local
	type (record), pointer :: CurrentRecord
	type (field), pointer :: CurrentField

	CurrentRecord => ThisRecord
	do while(associated(CurrentRecord))
!	  write(*,*)DescText, ' WriteRecord - a new record '
	  CurrentField => CurrentRecord%firstfield
	  do while(associated(CurrentField))
	    call WriteField('Next field in record', CurrentField)
	    CurrentField => CurrentField%next 
	  end do
	  CurrentRecord => CurrentRecord%next
	end do	

	return
	end subroutine WriteRecord
	
	subroutine WriteField(DescText, ThisField)
! Input
	character(len=*), intent(in) :: DescText
	type (field), pointer :: ThisField
! Local
	character*256 OutputText
	integer LenOutput

10017   format(2x,3a)

	if (associated(ThisField)) then
	  write(*,10017) DescText, '%data: ', trim(ThisField%data)
	  if(associated(ThisField%ident)) then
	    write(*,10017) DescText, '%ident%name: ', trim(ThisField%ident%name)
	  end if
	  call WriteArguments(ThisField%arg, OutputText, LenOutput) 
	  write(*,10017) DescText, ' Arguments: ', OutputText(1:LenOutput)
	end if

	return
	end subroutine WriteField

	subroutine GPSetName2Index(CurrentSetName, CurrentIndex)
! [01/01/04] Extracted from routine GPSet2Index
! Input
	character (len=*), intent(in) :: CurrentSetName
! Output
	character (len=*), intent(out) :: CurrentIndex
! Local
	character (len=llindex) :: IndexRoot, IndexRoot1, AddString, TryIndex
	type(symbol), pointer :: CurrentSymbol
	logical :: IsASetName, Finished
	integer :: i
!

! Truncate IndexRoot to at most LLINDEX-3 characters
	IndexRoot1 = CurrentSetName
	IndexRoot = IndexRoot1(1:llindex-3)

! Add integers to the end until find something which is not the name of a set.
	i = 1
	Finished = .false.
	do while(i.le.99 .and. .not.Finished)
	  AddString = PosintToString(i)
	  TryIndex = trim(IndexRoot)//'_'//AddString
!??	  TryIndex = trim(IndexRoot)//AddString
! Is this the name of a set? Go through the list of symbols in idlist to check.

! Really should also check whether this is the name of a currently active index. **
! To do that would need to keep a list of those and keep the list up-to-date.

	  IsASetName = .false.
	  CurrentSymbol => idlist
	  do while(associated(CurrentSymbol).and. (.not.IsASetName))
	    if(CurrentSymbol%type.eq.'set') then
	      if(EqualStringsIgnoreCase(TryIndex,CurrentSymbol%name)) then
	        IsASetName = .true.
	      end if
	    end if
	    CurrentSymbol => CurrentSymbol%next
	  end do ! while(associated(CurrentSymbol))

	  if(.not.IsASetName) then
	    CurrentIndex = TryIndex
	    Finished=.true.
	  end if
	
! Increment I
	  i = i+1
	end do ! while(i.le.99 .and. .not.Finished)

	return
	end subroutine GPSetName2Index

	subroutine GPSet2Index(CurrentSet, CurrentIndex)
! Input
	type (symbol), pointer :: CurrentSet
! Output
	character (len=*), intent(out) :: CurrentIndex
! Local
	character (len=llindex) :: IndexRoot, IndexRoot1, AddString, TryIndex
	type(symbol), pointer :: CurrentSymbol
	logical :: IsASetName, Finished
	integer :: i
!

! The CurrentSet passed here may be a quoted element. If so, set CurrentIndex
!  to equal CurrentSet
	if(CurrentSet%name(1:1).eq.'"') then
	  CurrentIndex = CurrentSet%name
	  return
	end if

! Need to check that the index we output is not the name of a set. [If it
!   is, TABLO will object.]
! Also need to check that the name is not too long.

! [01/01/04] Call routine GPSetName2Index
	call GPSetName2Index(CurrentSet%name, CurrentIndex)

	return
	end subroutine GPSet2Index

! Next replaced by a test running through GPDeclarationList
! [15/03/04] Call this HasBeenDeclaredGPOld to distinguish from the newer one
!	logical function HasBeenDeclaredGPOld(ThisSymbol)
! Tells if ThisSymbol has already been declared
! Input
!	type (symbol), pointer :: ThisSymbol
! Local
!	type (symbol), pointer :: TestSymbol
!	character(len=31) ThisName, TestName

!	HasBeenDeclaredGPOld = .false.

!	if(.not.associated(ThisSymbol)) then
!	  write(*,*) 'HasBeenDeclaredGPOld - symbol is not associated'
!	  return
!	end if
!	ThisName = trim(ThisSymbol%name)
!	call lower(ThisName)

! Run through all the symbols (they are in idlist) 
!  to find this one and check if it is declared
!	TestSymbol => idlist ;
!	do while(associated(TestSymbol))
!	  TestName = trim(TestSymbol%name)
!	  call lower(TestName)
! Are the two names the same ?
!	  if (TestName.eq.ThisName) then
!	    if(TestSymbol%DeclaredGP) then
!	      HasBeenDeclaredGPOld = .true.
!	      write(*,*)'HasBeenDeclaredGPOld:',trim(ThisName),' has been declared'
!	      Exit
!	    else
!	      write(*,*)'HasBeenDeclaredGPOld:',trim(ThisName),' says not declared'
!	    end if
!	  end if
!	  TestSymbol => TestSymbol%next
!	end do 

!	return
!	end function HasBeenDeclaredGPOld

! [15/03/04] No longer need this routine
!	subroutine MarkAsDeclaredGP(ThisSymbol)
! Indicate that ThisSymbol has been declared
! Input and output
!	type (symbol), pointer :: ThisSymbol
! Local
!	type (symbol), pointer :: TestSymbol
!	character(len=31) ThisName, TestName

!	ThisName = trim(ThisSymbol%name)
!	call lower(ThisName)

! Run through all the symbols (they are in idlist) 
!  to find this one and then mark it as declared
!	TestSymbol => idlist ;
!	do while(associated(TestSymbol))
!	  TestName = trim(TestSymbol%name)
!	  call lower(TestName)
! Are the two names the same ?
!	  if (TestName.eq.ThisName) then
!	    TestSymbol%DeclaredGp = .true .
!	    write(*,*)'MarkAsDeclaredGP:',trim(ThisName),' now declared'
!	  end if
!	  TestSymbol => TestSymbol%next
!	end do 

! Temporary
! Run through all the symbols again 
!  to find this one and then check it is shown as declared
!	TestSymbol => idlist ;
!	do while(associated(TestSymbol))
!	  TestName = trim(TestSymbol%name)
!	  call lower(TestName)
! Are the two names the same ?
!	  if (TestName.eq.ThisName) then
!	    if(TestSymbol%DeclaredGp) then
!	      write(*,*)'MarkAsDeclaredGP:',trim(ThisName),' shown as declared'
!	    else
!	      write(*,*)'MarkAsDeclaredGP:',trim(ThisName),' NOT shown as declared'
!	    end if
!	  end if
!	  TestSymbol => TestSymbol%next
!	end do 

!	return
!	end subroutine MarkAsDeclaredGP

	subroutine ArgToAllsGP(TheseArgs, OutputString)
! Input
	type (domain), pointer :: TheseArgs
! Output
	character (len=*), intent(out) :: OutputString
! Local
	integer LenOutputString
	type (domain), pointer :: CurrentArgument
	type (symbol), pointer :: CurrentSet
	character (len=128) CurrentIndex, CurrentAll
	integer LengthCurrentAll

	OutputString = ''
	LenOutputString = 0 

	CurrentArgument => TheseArgs
! The ALLs
	do while(associated(CurrentArgument))
	  CurrentSet => CurrentArgument%s
	  call GPSet2Index(CurrentSet, CurrentIndex)
	  CurrentAll = '(All,'//trim(CurrentIndex)//','//trim(CurrentSet%name)//')'
	  LengthCurrentAll = 7 + LEN(trim(CurrentIndex)) + LEN(trim(CurrentSet%name)) 
	  OutputString = OutputString(1:LenOutputString) // CurrentAll
	  LenOutputString = LenOutputString + LEN(trim(CurrentAll))
	  CurrentArgument => CurrentArgument%next
	end do

	return
	end subroutine ArgToAllsGP

	subroutine ArgToAllsGPUseIndexStack(TheseArgs, DimTheseArgs, &
		OutputString, AllSets, AllIndexes)
! This is used to get the quantifiers for the market clearing equation.
! The indexes used are different from those given by Set2Index.
!  The sets indexes range over and the indexes used are returned in 
!    arrays AllSets (the sets) and AllIndexes (the indexes). [These arrays
!   are like the relevant things in the index stack in the code for the program TABLO.]
! Input
	type (domain), pointer :: TheseArgs
	integer, intent(in) :: DimTheseArgs
! Output
	character (len=*), intent(out) :: OutputString
	character (len=*), intent(out) :: AllSets(DimTheseArgs)
	character(len=*), intent(out) :: AllIndexes(DimTheseArgs)
! Local
	integer LenOutputString
	type (domain), pointer :: CurrentArgument
	type (symbol), pointer :: CurrentSet
	character (len=128) CurrentIndex, CurrentAll
	integer LengthCurrentAll, NumAllsDone

	OutputString = ''
	LenOutputString = 0
	NumAllsDone = 0 

	CurrentArgument => TheseArgs
! The ALLs
	do while(associated(CurrentArgument))
	  NumAllsDone = NumAllsDone + 1
	  CurrentSet => CurrentArgument%s
	  call GPSet2Index(CurrentSet, CurrentIndex)
! For the present, just "mess" this about a bit. Later will probably need
!  to get more sophisticated
	  CurrentIndex=trim(CurrentIndex)//'@'
	  CurrentAll = '(All,'//trim(CurrentIndex)//','//trim(CurrentSet%name)//')'
	  LengthCurrentAll = 7 + LEN(trim(CurrentIndex)) + LEN(trim(CurrentSet%name)) 
	  OutputString = OutputString(1:LenOutputString) // CurrentAll
	  LenOutputString = LenOutputString + LEN(trim(CurrentAll))
! Add to output arrays
	  AllSets(NumAllsDone) = CurrentSet%name
	  AllIndexes(NumAllsDone) = CurrentIndex
	  CurrentArgument => CurrentArgument%next
	end do

	return
	end subroutine ArgToAllsGPUseIndexStack

	subroutine CoeffVarRealToGP(ThisCoeffVarReal, OutputString)
! Writes coeff/var followed by its arguments
! [24/03/04] Add the possibility of a real here
! Input
	type (GPCoeffVarReal) :: ThisCoeffVarReal
! Output
	character (len=*), intent(out) :: OutputString
! Local
	type (domain), pointer :: TheseArgs
	character (len=512) :: ArgString

	if(ThisCoeffVarReal%type.eq.'r') then
	  OutputString = ThisCoeffVarReal%name
	else
! Coeff or Variable
	  OutputString = ThisCoeffVarReal%name
	  if(ThisCoeffVarReal%type.eq.'c') then
! Upper case for Coefficients
	    call upper(OutputString)
	  end if
	  TheseArgs => ThisCoeffVarReal%arg
	  call ArgToGP(TheseArgs, ArgString)
	  if(ArgString.ne.'') then
	    OutputString = trim(OutputString)//'('//trim(ArgString)//')'
	  end if
	end if
	OutputString = trim (OutputString)

	return
	end subroutine CoeffVarRealToGP
	
	subroutine CoeffVarRealToGPUseIndexStack( &
! Input
	  ThisCoeffVarReal, DimIndexStack, AllSets, AllIndexes, &
! Output
	  OutputString)
! Writes coeff/var followed by its arguments
! [24/03/04] Add the possibility of a real here
! Input
	type (GPCoeffVarReal) :: ThisCoeffVarReal
	integer, intent(in) :: DimIndexStack
	character (len=*), intent(out) :: AllSets(DimIndexStack)
	character(len=*), intent(out) :: AllIndexes(DimIndexStack)
! Output
	character (len=*), intent(out) :: OutputString
! Local
	type (domain), pointer :: TheseArgs
	character (len=512) :: ArgString

	if(ThisCoeffVarReal%type.eq.'r') then
	  OutputString = ThisCoeffVarReal%name
	else
! Coeff or Variable
	  OutputString = ThisCoeffVarReal%name
	  if(ThisCoeffVarReal%type.eq.'c') then
! Upper case for Coefficients
	    call upper(OutputString)
	  end if
	  TheseArgs => ThisCoeffVarReal%arg
! Here need to look at index stack to see if take index from that
	  call ArgToGPUseIndexStack(TheseArgs, DimIndexStack, AllSets, AllIndexes, ArgString)
	  if(ArgString.ne.'') then
	    OutputString = trim(OutputString)//'('//trim(ArgString)//')'
	  end if
	end if
	OutputString = trim (OutputString)

	return
	end subroutine CoeffVarRealToGPUseIndexStack

	subroutine CoeffVarToGP(ThisCoeffVar, OutputString)
! Writes coeff/var followed by its arguments
! Input
	type (GPCoeffVarReal) :: ThisCoeffVar
! Output
	character (len=*), intent(out) :: OutputString

! [24/03/04] Just call routine CoeffVarRealToGP
	call CoeffVarRealToGP(ThisCoeffVar, OutputString)

	return
	end subroutine CoeffVarToGP

	subroutine ArgToGP(TheseArgs, OutputString)
! Write TheseArgs in form suitable for GEMPACK. eg "s1,f1"
! No start "(" or end ")" is added here
! Input
	type (domain), pointer :: TheseArgs
! Output
	character (len=*), intent(out) :: OutputString
! Local
	integer NumArg
	type (domain), pointer :: CurrentArgument
	type (symbol), pointer :: CurrentSet
	character (len=128) CurrentIndex, CurrentAll

	OutputString = ''
	NumArg = 0 

	CurrentArgument => TheseArgs
	do while(associated(CurrentArgument))
	  if(NumArg.gt.0) then
! comma to separate arguments
	    OutputString = trim(OutputString)//','
	  end if
	  CurrentSet => CurrentArgument%s
	  call GPSet2Index(CurrentSet, CurrentIndex)
	  OutputString = trim(OutputString)//trim(CurrentIndex)
	  NumArg = NumArg + 1
	  CurrentArgument => CurrentArgument%next
	end do

	return
	end subroutine ArgToGP

	subroutine ArgToGPUseIndexStack( &
! Input
	  TheseArgs, DimIndexStack, AllSets, AllIndexes, &
! Output
	  OutputString)
! Write TheseArgs in form suitable for GEMPACK. eg "s1,f1"
! No start "(" or end ")" is added here
! Input
	type (domain), pointer :: TheseArgs
	integer, intent(in) :: DimIndexStack
	character (len=*), intent(in) :: AllSets(DimIndexStack)
	character(len=*), intent(in) :: AllIndexes(DimIndexStack)
! Output
	character (len=*), intent(out) :: OutputString
! Local
	integer NumArg
	type (domain), pointer :: CurrentArgument
	type (symbol), pointer :: CurrentSet
	character (len=128) CurrentIndex, CurrentAll
	logical :: SetInIndexStack
	integer :: i

	OutputString = ''
	NumArg = 0 

	CurrentArgument => TheseArgs
	do while(associated(CurrentArgument))
	  if(NumArg.gt.0) then
! comma to separate arguments
	    OutputString = trim(OutputString)//','
	  end if
	  CurrentSet => CurrentArgument%s

! Need to see if this set is on the index stack. If so, use the index there.
	  SetInIndexStack = .false.
	  do i=1,DimIndexStack
	    if(EqualStringsIgnoreCase(CurrentSet%name, AllSets(i))) then
	      SetInIndexStack = .true.
	      CurrentIndex = AllIndexes(i)
	    end if
	  end do
	  if(.not.SetInIndexStack) then
!		[As before]
	    call GPSet2Index(CurrentSet, CurrentIndex)
	  end if

	  OutputString = trim(OutputString)//trim(CurrentIndex)
	  NumArg = NumArg + 1
	  CurrentArgument => CurrentArgument%next
	end do

	return
	end subroutine ArgToGPUseIndexStack

	subroutine DeclareGPLevelsVar(ThisGPCoeff, ThisGPLinearVar)
! Declare a levels variable (ThisGPCoeff) with the name of the associated 
!  linear variable (ThisGPLinearVar) stated explicitly in the declaration
! Input
	type (GPCoeffVarReal), intent(in) :: ThisGPCoeff, ThisGPLinearVar
! Local
	type (field), pointer :: ThisField
	type (symbol), pointer :: ThisSymbol
	character (len=80) CoeffString, VarString, OutputString
	logical :: IsChangeVar, LinVarHasBeenDeclared

	IsChangeVar = (ThisGPLinearVar%VarType.eq.'c')
	OutputString = 'Variable(Levels,'
	if(IsChangeVar) then
	  OutputString = trim(OutputString)//'Change,'
	end if
! Has this linear variable already been declared?
	LinVaRHasBeenDeclared = HasBeenDeclaredGP('vc', ThisGPLinearVar%name)
	if(LinVarHasBeenDeclared) then
!  Use Linear_Var=
	  OutputString = trim(OutputString)//' Linear_Var='//trim(ThisGPLinearVar%name)//')'
	else
!  Use Linear_Name=
	  OutputString = trim(OutputString)//' Linear_Name='//trim(ThisGPLinearVar%name)//')'
	end if

! Add the rest

	call newfield(ThisField)
	call newsymbol(ThisSymbol)
	ThisSymbol%name = ThisGPCoeff%name
	ThisField%ident => ThisSymbol
	ThisField%arg => ThisGPCoeff%arg

! If are declaring a levels variable, it is an internal error if this has already been declared as a Coefficient
	if (HasBeenDeclaredGP('c',ThisGPCoeff%name)) then
	  call InternalErrorStop(trim(ThisGPCoeff%name)//' has already been declared as a Coefficient')
	end if

! Declare it
	call DeclareFieldGP(OutputString, ThisField, ThisGPCoeff%labelinfo)

! Indicate that the linear variable has been declared
! [For the coefficient, this is done in routine DeclareFiledGP]
	if(.not.LinVarHasBeenDeclared) then
	  call AddGPDeclaration('vc', trim(ThisGPLinearVar%name))
	end if

	return
	end subroutine DeclareGPLevelsVar


	subroutine DeclareGPCoeffVar(ThisGPCoeffVar)
! Input
	type (GPCoeffVarReal), intent(in) :: ThisGPCoeffVar
! Local
	type (field), pointer :: ThisField
	type (symbol), pointer :: ThisSymbol
	character (len=80) CoeffVarString

	call newfield(ThisField)
	call newsymbol(ThisSymbol)
	ThisSymbol%name = ThisGPCoeffVar%name
	ThisField%ident => ThisSymbol
	ThisField%arg => ThisGPCoeffVar%arg
	if(ThisGPCoeffVar%type.eq.'c') then
	  if(ThisGPCoeffVar%IsParam) then	  
	    CoeffVarString = 'Coefficient (Parameter)'
	  else
	    CoeffVarString = 'Coefficient'
	  end if
	else if(ThisGPCoeffVar%type.eq.'v') then
	  if(ThisGPCoeffVar%OrigLevel.ne.' ') then
	    if(ThisGPCoeffVar%VarType.eq.'c') then
	      CoeffVarString = 'Variable (Change,Orig_Level='//trim(ThisGPCoeffVar%OrigLevel)//')'
	    else
	      CoeffVarString = 'Variable (Orig_Level='//trim(ThisGPCoeffVar%OrigLevel)//')'
	    end if
	  else
	    if(ThisGPCoeffVar%VarType.eq.'c') then
	      CoeffVarString = 'Variable (Change)'
	    else
	      CoeffVarString = 'Variable'
	    end if
	  end if
	else
	  call InternalErrorStop( &
	    'Unexpected value of ThisGPCoeffVar in routine DeclareGPCoeffVar. Value is: '//ThisGPCoeffVar%type)
	end if

	call DeclareFieldGP(CoeffVarString, ThisField, ThisGPCoeffVar%labelinfo)
 
	return
	end subroutine DeclareGPCoeffVar	

	subroutine DeclareFieldGP(DeclarationStart, ThisField, ThisLabellingInfo)
! Input
	character (len=*), intent(in) :: DeclarationStart  ! 'variable', 'coefficient' etc
!		DeclarationStart may contain qualifiers, eg (Parameter)
	type (field), pointer :: ThisField
	character (len=*), intent(in) :: ThisLabellingInfo
! Local
	character(len=512) OutputString, OutputStringAlls
	character(len=512) OutputStringArgs, CoeffVarName, ModDeclarationStart
	integer :: LenStart, LenSoFar, LenAlls
	character (len=3) TypeStart
	character (len=2) ThisType
	logical IsCoeff

10005   format(a)
10017   format(2x,3a)

! [22/03/04] Are the next 2 lines about ThisField relevant now??
! If not associated, exit
	if(.not.associated(ThisField)) Return
! If identifier is not associated, exit
	if(.not.associated(ThisField%ident)) Return

! Use ModDeclarationStart in place of DeclarationStart since will modify
!  it when declaring a levels variable
	ModDeclarationStart = DeclarationStart

!  Set up IsCoeff and ThisType
	TypeStart = DeclarationStart
	call lower(TypeStart)
! If this is declaring a levels variable, write out DeclarationStart and
!  then proceed as if declaring a coefficient
	if(TypeStart.eq.'var' .and. index(DeclarationStart, '(Levels').gt.0) then
	  call WriteOneLineToTAB(DeclarationStart)
! Reset ModDeclarationStart in this case
	  ModDeclarationStart = ''
! Now proceed as if declaring a coefficient
	  TypeStart = 'coe'
	end if
	  
	IsCoeff = .false.
	if(TypeStart.eq.'coe') then
	  IsCoeff = .true.
	else if(TypeStart.eq.'var') then
	end if
	if (IsCoeff) then
	  ThisType = 'cs'
	else
	  ThisType = 'vc'
	end if

! If already declared, exit
! WAS	if (HasBeenDeclaredGPOld(ThisField%ident)) Return
	if (HasBeenDeclaredGP(ThisType,ThisField%ident%name)) Return

	LenStart = LEN_TRIM(ModDeclarationStart)

	OutputString = ''

	OutputString=trim(ModDeclarationStart)//' '
	LenSoFar = LenStart+1

	call ArgToAllsGP(ThisField%arg, OutputStringAlls)
	OutputStringAlls = trim(OutputStringAlls)
	LenAlls = LEN_TRIM(OutputStringAlls)
	if(OutputStringAlls.ne.'') then
! Output Start now if ncessary
	  if(LenStart.gt.25) then
	    call WriteOneLineToTAB(OutputString)
	    OutputString = ' '
	  end if
	  OutputString = trim(OutputString)//' '//OutputStringAlls
	end if
	LenSoFar = LEN_TRIM(OutputString)

! Output a line now if ncessary
	if(LenSoFar.gt.40) then
	  call WriteOneLineToTAB(OutputString)
	  OutputString = ' '
	end if

! The name
	CoeffVarName = ThisField%ident%name
! Coefficients in upper case
	if(IsCoeff) then
	  call upper(CoeffVarName)
	end if
!	write(*,*)'ThisFile%ident%name): ', ThisField%ident%name
	OutputString = trim(OutputString)//' '//trim(CoeffVarName)
! The arguments
	if(OutputStringAlls.ne.'') then
! Opening (
	  OutputString = trim(OutputString)//'('
	end if
	call ArgToGP(ThisField%arg, OutputStringArgs)
	OutputString = trim(OutputString)//trim(OutputStringArgs)
	if(OutputStringAlls.ne.'') then
! Closing (
	  OutputString = trim(OutputString)//')'
	end if

! Labelling info and semicolon - on separate line - first output what we have
	call WriteOneLineToTAB(OutputString)
! Write the declaration
!	write(*,10017) 'Declaration: ', trim(OutputString)
! Labelling info and semicolon
	call PrepareLabelInfoForTABFile(ThisLabellingInfo, OutputString)
	OutputString = trim(OutputString)//' ;'
! Write to TAB file
	call WriteOneLineToTAB(OutputString)
! Mark this symbol as declared
! [15/03/04] No longer need this
!	call MarkAsDeclaredGP(ThisField%ident)
! Mark it in GPDeclarations
	call AddGPDeclaration(ThisType, trim(ThisField%ident%name))
 		
	return
	end subroutine DeclareFieldGP


	subroutine WriteNest(DescText, ThisNest)
	character*(*) DescText
	type (nest), pointer :: ThisNest
! Local
	type (netputlist), pointer :: elements
	type (netput), pointer :: ThisNetPut
	type (field), pointer :: IOField, qref, pref

10015   format(1x,a)

	elements => ThisNest%elements
	do while (associated(elements)) 
	  ThisNetPut => elements%np
	  IOField => ThisNetPut%f   ! I: or O: field
!	  if (ThisNetPut%output) then
!	    write(*,10015) ' An output'
!	  else
!	    write(*,10015) ' An input'
!	  end if
	  qref => ThisNetPut%qref
	  pref => ThisNetPut%pref

	  call WriteField('IOField', IOField)
	  if(associated(IOField)) then
!	    call DeclareFieldGP('Variable', IOField)
	  end if

	  call WriteField('QREF', qref)
	  if(associated(qref)) then
!	    call DeclareFieldGP('Coefficient', qref)
	  end if
	  call WriteField('PREF', pref)
	  if(associated(pref)) then
!	    call DeclareFieldGP('Coefficient', pref)
	  end if
	  elements => elements%next
	end do
	
	return
	end subroutine WriteNest

!!	Allocate a new GPDeclaration and assign default values
!!	for all the characteristics

	subroutine NewGPDeclaration(ThisGPDeclaration)
	type (GPDeclaration), pointer :: ThisGPDeclaration
	nullify(ThisGPDeclaration)
	allocate(ThisGPDeclaration)
	ThisGPDeclaration%type = ''  ! impossible type
	ThisGPDeclaration%name = ''  ! impossible name
	nullify(ThisGPDeclaration%next)
	return
	end subroutine NewGPDeclaration

!! Introduce a new GPDeclaration in the master list.

        subroutine AddGPDeclaration(type, name)      

        character (len=*), intent(in) :: type
        character (len=*), intent(in) :: name

        type (GPDeclaration), pointer :: ThisGPDeclaration

!	write(*,*)'AddGPDeclaration: type=', type

        if (.not.associated(GPDeclarationList)) then
!	  write(*,*)'AddGPDeclaration: not associated'
	  call NewGPDeclaration(GPDeclarationList)
          ThisGPDeclaration => GPDeclarationList
        else
!	  write(*,*)'AddGPDeclaration: is associated'
	  call NewGPDeclaration(GPDeclarationLast%next)
          ThisGPDeclaration => GPDeclarationLast%next
        end if
        ThisGPDeclaration%type = type
        ThisGPDeclaration%name = name
        nullify(ThisGPDeclaration%next)
        GPDeclarationLast => ThisGPDeclaration
        return
        end subroutine AddGPDeclaration

	logical function HasBeenDeclaredGP(type, name)
! Tells if this has already been declared
! Input
        character (len=*), intent(in) :: type
        character (len=*), intent(in) :: name
! Local
	type (GPDeclaration), pointer :: CurrentGPDeclaration
	character(len=2) CurrentType, TestType
	character(len=LengthGPDeclarationName) CurrentName, TestName

	HasBeenDeclaredGP = .false.

	TestType = type
	call lower(TestType)
	TestName = trim(name)
	call lower(TestName)

! Run throught the list of declarations seeking a match
	CurrentGPDeclaration => GPDeclarationList
	do while(associated(CurrentGPDeclaration))
	  CurrentType = CurrentGPDeclaration%type
!	  write(*,*)'HasBeenDeclaredGP: CurrentType=', CurrentType
	  call lower(CurrentType)
	  CurrentName = CurrentGPDeclaration%name
	  CurrentName = trim(CurrentName)
!	  write(*,*)'HasBeenDeclaredGP: CurrentName=', CurrentName
	  call lower(CurrentName)
	  if(CurrentName.eq.TestName) then
	    if(CurrentType.eq.TestType) then
	      HasBeenDeclaredGP = .true.
	      return
	    end if
	  end if
	  CurrentGPDeclaration => CurrentGPDeclaration%next
	end do	
	
	return
	end function HasBeenDeclaredGP

	subroutine WriteExogenousListOnCMF

	type (GPName), pointer :: CurrentGPName
	character (len=512) :: OutputString

	write(iocmfgp, '(a)') ' '
	write(iocmfgp, '(a)') '! Closure '
	write(iocmfgp, '(a)') '! ------- '
	write(iocmfgp, '(a)') '! '

! Go through GPNameList and build up list of all that are exogenous
	OutputString = 'exogenous'
	CurrentGPName => GPNameList
	do while(associated(CurrentGPName))
	  if(CurrentGPName%type.eq.'v') then
	    if(CurrentGPName%ExogEndogStatus.eq.'x') then
! Add this to the list
	      if(LEN_TRIM(OutputString).gt.60) then
		  write(iocmfgp,'(a)') trim(OutputString)
	          OutputString = ' '
	      end if
	      OutputString = trim(OutputString)//' '//CurrentGPName%name
	    end if
	  end if

	  CurrentGPName => CurrentGPName%next
	end do ! while(associated(CurrentGPName))

! Complete the list
	OutputString = trim(OutputString)//' ;'
	write(iocmfgp,'(a)') trim(OutputString)

	write(iocmfgp,'(a)') 'Rest endogenous ;'

	return
	end subroutine WriteExogenousListOnCMF

	subroutine WriteCMFFile

! Local
	character (len=512) :: OutputString, TempString

	write(iocmfgp, '(a)') 'Auxiliary files = '//trim(tabfileroot)//' ;'
	write(iocmfgp, '(a)') 'File Input = '//trim(tabfileroot)//'.har ;'

! Use "<CMF>" for updated file name
! WAS	write(iocmfgp, '(a)') 'Updated File Input = '//trim(tabfileroot)//'.upd ;'
	write(iocmfgp, '(a)') 'Updated File Input = <cmf>.upd ;'

	write(iocmfgp, '(a)') 'Extrapolation Accuracy File = yes ;'
	write(iocmfgp, '(a)') 'Log File = yes ;'
	write(iocmfgp, '(a)') 'Method = Euler ;  ! or Gragg'
	write(iocmfgp, '(a)') 'Steps = 4 6 8 ;  ! change as needed'

! Write the Closure to CMF file
	call WriteExogenousListOnCMF

	write(iocmfgp, '(a)') ' '
	write(iocmfgp, '(a)') '! Shocks '
	write(iocmfgp, '(a)') '! ------ '
	write(iocmfgp, '(a)') '! Please add them here '

	write(iocmfgp, '(a)') ' '
	write(iocmfgp, '(a)') '! Verbal Description '
	write(iocmfgp, '(a)') '! ------------------ '
	write(iocmfgp, '(a)') 'Verbal Description = '
	write(iocmfgp, '(a)') '! [Please add it here (can be several lines)]'
	write(iocmfgp, '(a)') '  ;  ! End of Verbal Description'

	end subroutine WriteCMFFile

	subroutine DoConstraints
! Run through the constraints, adding the appropriate equations to TAB file

! Local
	type(function), pointer :: CurrentFunction
	type(functionGP), pointer :: CurrentFunction2GP
	type (endowmentGP), pointer :: CurrentEndowment2GP
	type(field), pointer :: FirstField, SecondField
	integer :: NumFuncDone, NumELinesDone, NumConstraintsDone
	character (len=128) :: ConstraintVariableName, ConstraintAbbreviation
	character (len=128) :: EquationStart, OutputString
	type (GPCoeffVarReal) :: CurrentRFieldCoeff, ConstraintRFieldVar, ConstraintRFieldCoeff
	logical :: RelevantELineFound, IsLinear, IsLevels

	NumFuncDone = 0
	NumConstraintsDone = 0
	CurrentFunction => functionlist
	do while(associated(CurrentFunction))
	  FirstField => CurrentFunction%firstrecord%firstfield
	  NumFuncDone = NumFuncDone + 1
	  if(EqualStringsIgnoreCase(FirstField%label,'$constraint')) then
	    write(*,*)'Found $constraint:'//trim(FirstField%data)
	    ConstraintVariableName = FirstField%identtext
	    write(*,*)'ConstraintVariabelName is:'//trim(ConstraintVariableName)
	    ConstraintAbbreviation = trim(FirstField%label)//':'//trim(FirstField%data)

!  Re-initialize 'pre' messages again at start of each function
	    call InitializeMessages('pre',0)
	    Messages%PreLine(1) = 'Processing "'//trim(ConstraintAbbreviation)//'"'
	    write(*,*) trim(Messages%PreLine(1))

! Need to find this variable in an r: field in an e: line
	    RelevantELineFound = .false.
	    CurrentFunction2GP => FunctionGPList
	    do while(associated(CurrentFunction2GP))

	      if(CurrentFunction2GP.IsDemandFunctionGP) then
		NumELinesDone = 0
		CurrentEndowment2GP => CurrentFunction2GP%FirstEndowmentGP
		do while(NumELinesDone.lt.CurrentFunction2GP%NumEFields.and. .not.RelevantELineFound)
		  if(CurrentEndowment2GP%HasRField) then
		    CurrentRFieldCoeff = CurrentEndowment2GP%EndowRFieldCoeff
		    if(EqualStringsIgnoreCase(CurrentRFieldCoeff%name,ConstraintVariableName)) then
		      ConstraintRFieldVar = CurrentEndowment2GP%EndowRFieldVar
		      ConstraintRFieldCoeff = CurrentEndowment2GP%EndowRFieldCoeff
		      RelevantELineFound = .true.
		    end if
		  end if ! if(CurrentEndowment2GP%HasRField) then
		  NumELinesDone = NumELinesDone + 1
		  CurrentEndowment2GP => CurrentEndowment2GP%next
		end do ! while(NumELinesDone.lt.CurrentFunction2GP$NumEFields)
	      end if ! if(CurrentFunction2GP.IsDemandFunction) then

	      CurrentFunction2GP => CurrentFunction2GP%next
	    end do ! while(associated(CurrentFunction2GP))

	    if(.not.RelevantElineFound) then
	      call InternalErrorStop('Did not find variable after "'//trim(ConstraintAbbreviation)//'" in an r: field')
	    end if ! if(.not.RelevantElineFound) then

	    NumConstraintsDone = NumConstraintsDone + 1

	    if(NumConstraintsDone.eq.1) then
! Heading
	      call WriteOneLineToTAB(' ')
	      call WriteOneLineToTAB('![[! The $constraint Equations !]]!')
	    end if

! Echo the $constraint line to TAB file
	    call WriteOneLineToTAB(' ')
	    OutputString = '! $constraint:'//trim(FirstField%data)//' !'	    
	    call WriteOneLineToTAB(OutputString)
! Add the constraint equation to the TAB file
	    SecondField => CurrentFunction%firstrecord%next%firstfield
	    write(*,*)'SecondField%label is:', trim(SecondField%label)
	    write(*,*)'SecondField%data is:', trim(SecondField%data)
! Expect SecondField%data to begin "(linear)" or "(levels)"
	    IsLinear = EqualStringsIgnoreCase(SecondField%data(1:8), '(linear)')	
	    IsLevels = EqualStringsIgnoreCase(SecondField%data(1:8), '(levels)')
	    if(IsLinear) then
		EquationStart = 'Equation (Linear)'
	    else if(IsLevels) then
		EquationStart = 'Equation (Levels)'
	    else
		call ErrorStop('Expect "(Linear)" or "(Levels)" at start of $constraint equation')
	    end if
	    OutputString = trim(EquationStart)//' E_'//trim(ConstraintVariableName)
! Write this line
	    call WriteOneLineToTAB(OutputString)
! Write the rest of the equation exactly as in the MGE file - skip (linear) or (levels)
	    OutputString = SecondField%data(9:)
	    call WriteOneLineToTAB(OutputString)
! Mark the auxiliary variable as endogenous
	    call MarkAsEndogenous(ConstraintRFieldVar)

	  end if ! if(EqualStringsIgnoreCase(FirstField%label,'$constraint')) then

	  CurrentFunction => CurrentFunction%next
	end do  ! while(associated(CurrentFunction))

!  Re-initialize 'pre' messages again at start of each function
	call InitializeMessages('pre',0)

	return
	end subroutine DoConstraints

	subroutine OmitAllOfSpecifiedString(SpecString, InputOutputString)
! Input
	character (len=*), intent(in) :: SpecString
! Input and output
	character (len=*) :: InputOutputString
! Local
	integer CurrentLength, n1, LenStringOmitting, end1pos, st2pos

	LenStringOmitting = LEN_TRIM(SpecString)
12	continue
	CurrentLength = LEN_TRIM(InputOutputString)
	n1 = INDEX(InputOutputString, SpecString)
	if(n1.gt.0) then
! END1POS is end of first string left after make this omission
	  end1pos=n1-1
! ST2POS is start of second string left after make this omission
	  st2pos=n1+LenStringOmitting
	  if(end1pos.ge.1.and.st2pos.le.CurrentLength) then
	    InputOutputString = InputOutputString(1:end1pos)//InputOutputString(st2pos:)
	  else if(end1pos.le.0.and.st2pos.le.CurrentLength) then
	    InputOutputString = InputOutputString(st2pos:)
	  else if(end1pos.ge.1.and.st2pos.gt.CurrentLength) then
	    InputOutputString = InputOutputString(1:end1pos)
	  else
	    InputOutputString = ''
	  end if
! Do it again
	  goto 12
	end if
	 
	return
	end subroutine OmitAllOfSpecifiedString

end module mpsGP